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Abstract 

A  computer  code  has  been  written  which  calculates  the  small  deformation  stress 
and  strain  fields  of  a  medium  consisting  of  a  pack  of  rigid  spherical  particles  embedded 
in  an  elastic  Hookean  matrix.  The  stress  and  strain  tensors  can  be  calculated  at  any  point 
in  the  medium  to  within  a  user-specified  accuracy.  Average  mechanical  properties  of  the 
medium  are  also  output  by  the  code.  The  code  has  been  used  to  simulate  systems 
consisting  of  thousands  of  particles  in  a  finite  pack.  Optionally,  the  code  treats  an  infinite 
pack  made  up  of  repeating  3D  rectangular  cells  of  a  particle  pack.  The  multipole 
expansion  technique  used  to  solve  the  equations  of  small  defonnation  for  the  elastic 
medium  consists  of  truncated  sums  of  complete  orthogonal  vector  spherical  harmonics. 
Techniques  are  presented  which  improve  the  convergence  of  the  solution  when  the 
particles  are  in  close  proximity  for  highly  filled  materials.  The  code  has  been  tested 
against  exact  solutions  of  configurations  consisting  of  a  few  particles  as  well  as  infinite 
packs  of  particles  in  body-centered  cubic,  face-centered  cubic,  and  simple  cubic  lattice 
arrangements.  The  code  has  been  used  to  estimate  mechanical  properties  of  a  variety  of 
monomodal  and  bimodal  particle  packs  of  different  packing  densities  for  a  matrix 
material  with  a  variety  of  elastic  constants. 
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1.  Introduction 


Mechanical  models  of  particulate  materials  can  be  categorized  in  several  ways. 

•  Particle  properties:  isotropic,  rigid  particles,  deformable  particles,  voids 

•  Particle  shape  (Davis,  1999):  spherical  (Goodier,  1933;  Walpole,  1972),  ellipsoidal 
(Robinson,  1951;  Edwards,  1951;  Eshelby,  1957;  Buchalter,  1994),  irregular 
(Yamamoto,  et  al,  1999),  chopped  fiber  (Lopez -Pamies,  1996) 

•  Matrix  properties:  isotropic,  linearly  elastic,  nonlinear  elastic  (Chen,  et  al.,  1993), 
nonlinear  viscoelastic  (Schapery,  1986)  ;Matous,  et  al.,  2005),  nonlinear  viscoplastic 
(Olsen  et  al.,  1998) 

•  Particle-matrix  bond:  bond  strength  less  than  (Gent,  1980)  or  greater  than  (Gent  et 
al.,  1984)  matrix  failure  strength  (T.  Smith,  1959;  Sangani,  1997) 

•  Microstructure  evolution:  microstructure  changes  (e.g.,  particle  debonding)  or  does 
not  change  with  application  of  boundary  condition  history  such  as  loads, 
temperatures,  and  chemical  environments 

•  Computational  efficiency:  empirical  models  (van  der  Poel,  1958;  Hashin,  1964;  J. 
Smith,  1974,  1975,  1976),  semi-empirical  and  bounded  models  (Hashin,  1983; 
Christensen,  1990;  Aravas,  2005),  special  case  models  (Sangani,  1987),  physics- 
based  models  with  simplified  microstructural  input  (Christensen,  1990;  Vratsanos, 
1993),  intensive  physics-based  models  with  detailed  microstructural  input  (Davis,  et 
al.,  1993;  Phan-Thien  et  al.,  1994;  Torquato,  1997) 

Many  users  of  particulate  material  models  need  mechanical  constitutive  models 

that  can  be  placed  into  a  finite  element  code  to  predict  the  integrity  of  structures  partially 

or  wholly  comprised  of  particulate  composites.  The  constitutive  relation  should  satisfy 
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three  criteria:  (1)  reasonable  accuracy,  (2)  computational  tractability,  and  (3)  numerical 
stability.  Unfortunately,  computational  tractability  and  reasonable  accuracy  are  difficult 
to  simultaneously  obtain.  The  more  microstructural  information  one  puts  into  a  model, 
the  more  accurate  it  can  be,  but  the  more  computationally  intensive  it  will  become. 

The  construction  of  a  routinely  usable  particulate  model  therefore  becomes  a  search 
for  those  physical  phenomena  that  have  greatest  impact  for  the  least  computational 
burden.  In  order  to  filter  through  the  most  influential  physical  phenomena,  we  must  have 
models  that  may  themselves  be  too  detailed  and  computationally  intensive  for  practical 
use,  but  can  guide  and  validate  simpler  models.  With  a  detailed  model,  one  can  quantify 
the  accuracy  and  computational  or  stability  gain  for  each  simplifying  assumption.  Broad 
treatments  with  many  references  on  particulate  theories  and  modeling  issues  may  be 
found  in  Christensen  (1979),  Cusak  (1987),  and  Torquato  (2002). 

In  this  paper,  we  describe  a  computationally  intensive  model  for  rigid,  bonded, 
spherical  particles  in  an  isotropic,  linearly  elastic  matrix  using  Navier  multipoles.  Later 
papers  will  present  our  efforts  at  computationally  tractable  statistical  mechanical  models 
condensed  from  the  detailed  treatment  of  this  paper  and  with  added  microstructural  detail 
such  as  debonding  particles,  viscoelastic  and  large  deformation  matrix  materials,  and 
nonspherical  particles. 

When  the  authors  first  employed  a  multipole  expansion  technique  to  treat  the 
mechanical  response  of  a  particulate  material  (Davis,  et  ah,  1993),  several  dozen  particles 
could  be  efficiently  modeled  with  Navier  multipoles.  Since  that  early  effort,  another 
project  was  launched  to  refine  and  enhance  the  capabilities  of  the  method.  A  more 
elegant  mathematical  foundation  based  on  vector  spherical  harmonics  was  adopted  that 
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greatly  simplified  the  multipole  expansion  method.  This  formalism  led  to  the 
straightforward  implementation  of  addition  theorems  that  give  the  solutions  in  terms  of  a 
large  number  of  linear  algebraic  equations.  The  use  of  addition  theorems  allows  several 
efficiency  measures  to  be  incorporated  into  the  solution  of  the  problem  so  that  thousands 
of  particles  in  a  highly  filled  medium  can  be  efficiently  modeled  on  current  computers. 
An  efficient  technique  was  devised  to  treat  the  interaction  between  particles  of  very  close 
proximity  with  very  high  order  multipole  tenns. 

The  original  paper  treated  a  finite  collection  of  particles.  The  method  has  since 
been  extended  to  approximate  an  infinite  medium  with  repeating  3D  rectangular  cells  of 
particles.  The  mathematical  methods  described  in  this  paper  are  programmed  in  the 
PARMECH  computer  code. 

Independently,  work  was  carried  along  similar  lines  by  other  researchers.  Kushch 
(1997)  used  a  multipole  expansion  technique  to  model  aligned  spheroidal  particles  (not 
necessarily  rigid)  in  periodic  repeating  unit  cells.  Addition  theorems  were  used  to  reduce 
the  problem  to  an  algebraic  system  of  equations.  Numerical  solutions  were  given  for  a 
single  particle  in  a  unit  cell  corresponding  to  a  periodic  lattice. 

Sangani  and  Mo  (1997)  also  used  a  multipole  expansion  technique  to  treat 
repeating  cells  of  spherical  particles.  Accurate  calculations  were  made  of  highly-filled 
systems  containing  up  to  32  random  particles  of  a  single  size  in  the  unit  cell. 

Despite  the  previous  work  in  this  area,  there  is  a  lack  of  calculations  of  model 
systems  that  can  be  used  to  guide  the  construction  of  particulate  constitutive  models.  In 
particular,  the  authors  found  no  systematic  calculations  treating  systems  where  the 
particles  differed  in  size  by  a  significant  amount.  Of  course  such  systems  are  challenging 
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because  of  the  large  number  of  small  particles  that  must  be  included  when  a  statistically 
significant  number  of  larger  particles  are  treated.  Nevertheless,  such  calculations  are 
critical  in  validating  models  that  treat  collections  of  particles  of  different  sizes. 

Most  of  the  calculations  presented  in  this  paper  are  for  a  matrix  material  with 
Poisson’s  ratio  near  Vi.  When  the  particles  are  in  close  proximity,  the  stress 
concentrations  between  close  particles  are  particularly  large  for  this  condition.  To  treat 
these  concentrations  accurately  with  a  multipole  expansion  technique  is  a  challenge  and 
highlights  the  efficiency  of  the  nearest  neighbor  treatment  that  allows  the  effective 
treatments  of  particles  in  close  proximity. 

In  this  paper  the  mathematical  fonnulation  based  on  vector  spherical  harmonics  of 
the  multipole  expansion  technique  is  summarized.  It  is  shown  how  the  addition  theorems 
can  be  applied  to  the  governing  equations  to  satisfy  the  appropriate  boundary  conditions. 
Efficiency  measures  implemented  into  the  code  are  explained  and  the  nearest  neighbor 
treatment  is  derived.  Next,  a  comparison  of  predictions  of  the  PARMECH  code  is  made 
with  exact  solutions  for  a  simple  two-particle  system  as  well  as  infinite  repeating  lattice 
arrangements  of  particles.  Predictions  of  the  code  for  random  arrangements  of  particles 
are  also  presented.  We  treat  in  detail  systems  containing  a  single  particle  size  and  also 
two  distinct  particle  sizes. 

2.  Elements  of  particulate  models 

Three  questions  that  should  be  addressed  when  building  a  microstructural 
particulate  model  are 

(1)  What  are  the  relevant  properties  of  the  individual  phases? 
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(2)  How  are  the  phases  geometrically  arranged  relative  to  each  other? 

(3)  What  are  the  laws  of  interaction  among  the  phases  when  subject  to  boundary 
conditions  of  interest? 

Within  the  scope  of  our  work,  we  will  answer  the  first  question  in  this  section,  the  second 
question  in  Section  3,  and  the  third  question  in  Section  4. 

In  the  particulate  composites  of  interest  to  us,  the  particle  Young’s  modulus  is 
orders  of  magnitude  higher  than  matrix  modulus.  Hence  we  treat  the  particles  as  rigid. 
However,  we  emphasize  that  our  multipole  formalism  works  equally  well  for  a  set  of 
nonrigid  particles,  each  with  its  own  Lame  constants.  The  algebra  involved  roughly 
doubles  but  the  conceptual  framework  is  no  more  complicated. 

We  treat  the  matrix  as  an  isotropic,  Hookean  solid  described  by  the  two  elastic 
Lame  constants,  X  and  p.  These  two  constants  are  input  to  the  model. 

We  will  treat  the  case  where  the  particles  are  initially  bonded  and  stay  bonded  to 
the  matrix.  Again,  the  Navier  multipole  formalism  can  treat  particles  that  debond  from 
the  matrix  after  a  certain  threshold  criterion  is  met.  However,  since  there  are  a  variety  of 
boundary  conditions  that  can  be  applied  to  debonding  particles,  there  is  much 
experimental  and  theoretical  work  still  needed  to  write  the  simulation  code  for  debonding 
particulate  materials.  We  have  worked  through  some  of  the  algebra  on  several  debonding 
scenarios  but  these  treatments  are  not  programmed  or  mature  enough  to  present  here.  We 
hope  to  present  them  in  subsequent  work.  It  is  clear,  however,  that  the  computational 
burden  will  increase  substantially,  more  than  an  order  of  magnitude,  when  debonding  is 
pennitted.  It  is  also  clear  that  the  nature  of  debonding  and  the  subsequent  particle-matrix 
boundary  conditions  for  a  given  particulate  material  is  not  trivial  to  characterize. 
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3.  Random  particle  packing  with  arbitrarily  broad  particle  size  distributions 

To  make  detailed  microstructural  calculations  of  particulate  materials,  one  must 
have  the  detailed  microgeometry  of  the  particle  pack.  Hence  a  particle-packing  algorithm 
must  be  available  which  can  place  a  sufficiently  large  number  of  particles  to  create  a 
faithful  statistical  sampling  of  a  full  pack.  Many  packing  algorithms  exist  (Davis,  1999). 
Two  important  categories  of  particle  packing  models  are  ballistic  deposition  and 
molecular  dynamics.  In  ballistic  deposition,  the  particles  are  added  to  the  pack  one  at  a 
time.  Each  particle  rolls  about  on  the  pack  until  it  finds  a  stable  equilibrium  point  resting 
upon  other  particles  of  the  pack  and/or  upon  a  simulated  container  wall.  In  molecular 
dynamics  (as  the  name  implies),  the  particles  are  placed  in  a  container  and  allowed  to 
move  and  collide  with  each  other  and  the  container  wall  as  the  wall  slowly  closes  them 
into  a  close  pack  configuration.  The  ballistic  deposition  algorithms  are  computationally 
much  faster  but  do  not  pack  to  as  high  a  packing  density  as  the  molecular  dynamics 
methods.  A  third  category  of  packing  algorithms  grows  the  particles  from  randomly 
distributed  seeds  (Kochevets  et  ah,  2001)  which  then  move  to  avoid  overlap  as  growth 
causes  them  to  bump  into  each  other. 

The  importance  of  the  packing  algorithm  can  hardly  be  overemphasized  when 
dealing  with  particulate  composites  that  are  near  close  pack.  In  close  pack,  the  particles 
are  nearly  touching,  that  is,  the  mean  surface-to-surface  distance  is  much  less  that  the 
particle  radii.  The  mechanical  properties  near  close  pack  are  ruled  by  the  large  stress 
concentration  spikes  that  occur  between  nearly  touching,  rigid  particles  (Figure  1).  The 
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height  and  size  of  the  stress  concentration  spikes  are  ruled  by  the  distribution  of  surface- 
to-surface  distances,  being  inversely  proportional  to  the  surface-to-surface  distance  for 


Figure  1.  Stress  concentration  spike  (trace  of  the  stress  tensor)  plotted  for  the  region  around  two  rigid 
spheres  embedded  in  a  nearly  incompressible  matrix  under  uniaxial  tension  in  the  z  direction.  These  spikes 
dictate  mechanical  properties  in  highly  filled  particulate  materials. 


nearly  incompressible  matrix  material.  From  experience,  we  have  learned  that  even 
though  a  packing  algorithm  may  give  the  correct  volume  fraction  of  particles,  if  it  does 
not  mimic  the  actual  surface-to-surface  distance  distribution  of  the  real  material,  the 
predicted  mechanical  properties  can  be  in  error  by  a  few  hundred  percent.  A 
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considerable  burden  is  therefore  placed  on  the  researcher  to  understand  how  the  real  pack 
was  constructed  and  adapt  packing  algorithms  that  yield  the  same  surface-to-surface 
distance  distributions  with  good  fidelity. 

The  particle  packs  of  interest  to  us  have  broad  size  distributions.  Most  packing 
algorithms  cannot  treat  broad  distributions  because  of  the  simple  problem  of  numbers. 
Consider  a  bimodal  pack  where  the  two  particle  sizes  differ  by  a  factor  of  one  hundred. 
Suppose  the  coarse  and  fine  modes  are  roughly  equivalent  in  volume  quantity.  Then, 
since  the  particle  diameters  differ  by  10',  the  particle  volumes  differ  by  10  .  If  one  must 
place  several  thousand  coarse  particles  to  get  a  reasonable  statistical  sampling  of  the 
pack,  it  follows  one  must  place  106x(several  thousand)  =  several  billion  small  particles  to 
fill  the  interstitial  spaces  of  the  larger  particles. 

To  circumvent  this  problem,  we  devised  a  scheme  called  a  reduced-dimension 
packing  algorithm  (Davis  et  ah,  1990),  the  short-comings  of  which  have  been  recently 
corrected  (Webb  et  al.,  2006)  to  provide  a  robust  packing  code  PARPACK  for  broader 
size  distributions  than  are  currently  available  in  other  packing  codes  (Figure  2).  On 
current  desktop  computers,  PARPACK  can  create  simple  packs  with  a  few  million 
particles  (far  more  than  is  generally  needed  to  obtain  the  statistical  flavor  of  simple 
packs)  in  a  few  minutes  and  more  complex  packs  of  the  same  size  in  a  few  hours. 
PARPACK  can  also  generate  coordination  numbers  for  the  modes  of  the  pack  and  their 
radial  distribution  functions.  Generating  higher-order  statistical  functions  would  be  a 
simple  matter  but  has  not  been  done  at  this  point. 

We  use  PARPACK  to  generate  packs  used  by  the  PARMECH  code  described  in 
this  paper. 
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Figure  2.  Visualization  of  a  reduced  dimension  pack.  The  smaller  particles  are  constrained  to  lie  close  to 
the  pack  axis.  Then  the  statistics  of  the  inner  small  particle  columns  are  projected  into  the  outer  part  of  the 
pack  where  the  larger  particles  place  themselves  according  to  small  particle  statistics.  Hence,  large 
particles  lying  in  the  outer  regions  of  the  pack  respond  to  small  particle  statistics  rather  than  billions  of 
small  particles  themselves.  In  this  way  we  mimic  full  3D  pack  statistics  with  broad  size  distribution  with  a 
tiny  fraction  of  the  particles  needed  in  a  3D  pack.  This  cut-away  picture  peers  past  some  large  particles 
into  the  small-particle  core  of  the  pack. 


4,  The  laws  of  interaction 

In  this  section  we  describe  the  mathematical  foundations  of  the  multipole  method 
of  solution  used  by  the  PARMECH  code.  In  subsequent  sections,  we  discuss  efficiency 
measures  employed  by  the  code  to  speed  the  solution  process.  We  also  summarize  the 
capabilities  of  the  code. 


10 


4. 1  General  solution  of  the  Navier  equation 

The  Navier  equation  for  elastic,  isotropic,  Hookean  material  is  given  by 
(A  +  2//)V(V  -u)  -  n  V  x  V  x  u  =  0  (1) 

where  A  and  //  are  the  Lame  constants.  The  small  deformation  displacement  vector  is 
denoted  by  u  .  Cartesian  components  will  be  used  for  all  displacement  vectors  so  that 
they  may  be  readily  added  when  the  interactions  of  many  spheres  are  included.  However, 
each  Cartesian  component  of  u  is  in  spherical  polar  coordinates  ( r,0,(p )  centered  on  the  i- 
th  sphere,  which  that  particular  displacement  vector  field  is  describing: 

u'k  =  u'k (r,  6,  (p)  ,  k  =  x,  y, z  .  (2) 

A  general  solution  to  the  Navier  equation  for  the  region  outside  a  sphere  i  of 


radius  a  is 


oo  L+ 1 


'  (F)  =  Z  X  X  Bl-nJuun  (7) 


(3a) 


L= 0 n=L- 1  m=-n 


Ui.nm  (r)  =  VL  nm  (?)  +  dn  (r 2  -  a2  )V[V  •  VL  nm  (F)] 


(3b) 


Yi.nJr)  = 


f-1 

vC 


L+ 1 


YnLm(0,<p) 


(3c) 


d„  = 


A  +  ju 


2  [A  +  2/u  +  (A  +  3  /u)n\ 


(3d) 
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(r2  -a2)V[V  •  VL  nm(r)\  = 


(3e) 


-  (2 n  - \)sjn(n  + 1) 

< 

0  ,  Z  =  «  or  n  + 1 

where  the  constants  B‘L  nm  are  called  the  multipole  moments  for  sphere  i  and  r  is  the 
distance  from  the  center  of  the  sphere  (r  >  a)  .  The  Y]nm  (6,  <p)  are  vector  spherical 
hannonics  defined  by  (Varshalovich  et  ah,  1988,  p.  210) 

fL(0,‘P)  =  'LCZuYl,A8.'P)et  (4) 

m\k 

where  the  unit  vectors  are  defined  by 
e_ i  =  (ex  -  zey)/21/2 

e0  =  ez  (5) 

e+i  =  -(ex  +  zey)/ 2 1/2  , 

the  C™,u  are  Clebsch-Gordon  coefficients,  and  the  YLm  are  scalar  spherical  hannonics. 
These  vector  spherical  hannonics  are  orthonormal 

2?r  1 

I  d<pf  d(oosff)  YZ  ■  f.t  =  SLL.S,,S„m.  ,  (6) 

0  -1 

and  complete,  in  that  there  always  exists  a  unique  set  of  coefficients  AL  nm  in  which  any 
arbitrary  piecewise  continuous  vector  function  F(6,cp)  can  be  expanded: 

oo  L+ 1  n 

wp)=££  Yal.JL(0,<p)  •  (7) 

L= 0 n=L- 1  m=—n 

The  Navier  solution  (3)  about  a  sphere  simplifies  considerably  at  r  =  a: 


f  \n  {  \n+2 

a  a 


Y. 


n+ 1 


,  L  =  n  - 1 
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(8) 


u{r) 


r\=a 


co  L+ 1  n 


II  I B[,JL(e.v)- 


The  completeness  of  the  vector  spherical  harmonics  implies  that  any  displacement  field 
obeying  the  Navier  equation  about  a  sphere  can  be  described,  from  a  simple  translational 
displacement  of  a  bonded  sphere  to  a  complex  displacement  for  an  unbonded  sphere. 
Only  the  bonded  sphere  is  treated  in  this  paper.  The  orthonormality  of  the  vector 
spherical  harmonics  allows  the  boundary  condition  about  a  sphere  to  be  applied  without 
the  need  for  time-consuming  numerical  integrations. 

The  general  solution  to  the  Navier  equation  for  N  spherical  particles  can  be 
written  as 


N  In 

u (r)  =  Yj u' '  (r  -  rt )  +  £  J  Fnm rYxnm  (0,<p),  (9) 

i= 1  n= 0  m=—n 

where  ri  is  the  position  of  particle  i.  The  first  sum  in  (9)  is  a  sum  of  single  particle 

solutions  defined  in  (3)  and  goes  to  zero  as  r  — >  oo.  The  second  sum  represents  the 
applied  field  far  from  the  particles  with  Fnm  derived  from  the  applied  strains.  The 
specific  solution  to  a  given  problem  is  therefore  obtained  by  determining  the  coefficients 
B‘l  nm  of  (3)  that  satisfy  the  boundary  conditions  for  the  problem. 

4.2  Boundary  conditions  and  addition  theorems 

The  PARMECH  code  currently  treats  the  problem  for  which  the  particles  are 
bonded  to  the  matrix  material.  The  boundary  condition  at  the  surface  of  a  particle  is 
given  by  a  rigid  displacement  and  rotation  of  the  particle  surface  and  can  be  written  as 


u 


(10) 
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where  the  first  term  in  the  sum  gives  the  translation  of  the  surface  and  the  second  term 
gives  the  rotation  of  the  surface. 

Large  computation  efficiency  is  achieved  through  the  powerful  mathematical 
relationship  of  addition  theorems  that  expresses  a  vector  spherical  hannonic  about  a  point 
r.  in  tenns  of  a  sum  of  vector  spherical  harmonics  about  a  point  ri .  This  allows  the 

general  solution  (9)  evaluated  at  the  surface  of  sphere  i  to  be  rewritten  in  terms  of  a  sum 
of  vector  spherical  harmonics  with  coordinates  centered  about  the  point  r.  as 


u\r)  =  Y, 

\r\=n 


L.n.m 


B‘  ,V  V  gi 

^ L,nm  '  ^ L  ,n  m  ^ L,nm 


y=l  L'  ,n'm 


Yi(0,<p) 


1  L+ 1  n  2  ri 

+i  Z  Z  i  z^,w  fii,jYL  (e,<p) 

L= 0  n=L—\  m=—n  n'= 0  m'=—ri 


where 


(11) 


oo  L+ 1  n 

z  =zz  z 


L,n,m  L= 0  n=L- 1  m=—n 


(12) 


and  the  off-center  expansion  coefficients  ’lj  and  are  derived  from  the 

addition  theorems  of  the  spherical  harmonics  and  depend  only  on  the  geometry  of  the 
particle  pack.  Their  derivation  is  given  in  Appendix  A. 

The  boundary  conditions  are  now  easily  applied  by  setting  (10)  equal  to  (1 1)  and 
applying  the  orthonormality  condition  of  (6),  then  rearranging  terms  to  obtain  the 
following  set  of  linear  equations: 


B 


L,nm 


Ll^nl 


■I  T.b(w„c, 


L'  ,ritri  ,ij 
L.nm 


7=1  L',n'm' 


y  7^  rn'm'Jj 
/  i  /  j  rim'  J  L.nm 


n'= 0  m'=—ri 


(13) 
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Given  the  displacements  ylm  and  rotations  \j/'m  ,  (13)  can  be  solved  for  the  coefficients 
B‘l  nm ,  thus  completely  specifying  the  displacement  vector  field  and  providing  a  complete 
solution  to  the  problem. 

The  solution  process  is  easily  modified  if  the  rotations  and  displacements  are  not 
specified  but  that  the  particles  are  in  equilibrium  (i.e.,  forces  and  torques  on  each  particle 
are  zero)  is  specified.  For  this  case  it  can  be  shown  that 

B '  lm  =  0  and  B[  itn  =  0 ,  (force  and  torque  on  particle  i  are  zero).  (14) 

The  displacements  y'm  and  rotations  y/'m  can  then  be  determined  from  (13). 

4.3  Efficiency  measures  and  nearest-neighbors  treatment 

When  solving  (13),  the  sums  can  be  truncated  at  some  maximum  value  Lmax 
depending  on  the  desired  accuracy  of  the  solutions.  It  is  important  to  truncate  this  sum  at 
the  lowest  reasonable  value  since  the  number  of  terms  Nmp  in  a  given  particle’s  multipole 


expansion  grows  as 

Nmp  =  3(Lmax+lf .  (15) 

When  stress  concentrations  in  the  problem  are  large,  as  they  are  for  packs  where  the 
particles  are  in  close  proximity,  a  large  value  of  Lmax  is  required  for  good  accuracy.  The 
values  of  Lmax that  are  required  for  a  given  accuracy  are  discussed  in  later  sections. 

Another  efficiency  measure  is  possible  by  looking  at  the  form  for  the  off-center 
expansion  coefficients 


(16) 
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where  Ry  is  the  distance  between  spheres  i  and  j.  Thus,  the  off-center  expansion 
coefficients  become  small  when  the  particles  are  far  apart  for  large  values  of  L  and L' 
which  allows  an  earlier  truncation  of  the  sum. 

A  further  important  efficiency  measure,  which  we  refer  to  as  the  nearest- 
neighbors  treatment,  is  included  in  the  code  to  allow  an  accurate  treatment  of  particles 
when  they  are  in  close  proximity.  This  approximation  is  of  value  when  the  Poisson’s 
ratio  is  near  0.5. 

The  greatest  stress  concentration  in  a  particle  pack  is  in  the  region  between  two 
nearly  touching  particles.  The  stress  concentrations  are  particularly  high  when  the 
Poisson’s  ratio  of  the  matrix  material  is  near  0.5  and  there  is  an  axial  displacement  of  the 
two  particles.  The  asymptotic  load  transfer  between  the  two  spheres  in  that  case  is 
proportional  to  s  ,  where  s  is  a  dimensionless  distance  between  the  sphere  surfaces.  On 
the  other  hand,  when  the  Poisson’s  ratio  is  not  near  0.5  or  for  a  shear  displacement 
between  particles  when  Poisson’s  ratio  is  near  0.5,  the  asymptotic  load  transfer  (Phan- 
Thien  et  ah,  1994,  Ch.  4)  is  proportional  to  khV1). 

These  two-body  interactions  can  be  described  with  a  large  number  of  multipoles 
separately  from  the  remainder  of  the  pack.  For  some  calculations,  while  the  general  pack 
used  an  Lmax  of  12,  the  nearest-neighbors  treatment  used  an  Lmax  of  200  for  the  two-body 
interaction  only. 

We  let  NNi  be  the  set  of  nearest  neighbors  of  particle  i  and  write  the  multipole 
expansion  coefficients  of  particle  i  as 


^ '  L.nm  y  LO^ fl\ 


+  V'm^L\^n\  + 


Ik 


J  +  A' 

L,nm  ^L.nm 


jeNNi 


(17) 
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The  first  two  terms  on  the  right  hand  side  of  (17)  describe  the  displacement  and  rotation 
of  particle  i,  the  B‘/jnm  represents  the  contribution  due  to  the  tensile  interaction  of  particle 

i  with  particle  j  and  must  be  treated  with  high  order,  and  the  A\  nm  give  the  remaining 

contribution  to  the  multipole  coefficients  and  need  not  be  treated  with  high  order. 

If  we  substitute  (17)  into  (13)  for  the  nearest  neighbors  of  particle  i  we  obtain  the 
following  relationship: 


I 

jsNNi 


D'J  ,  A1 
UL,nm  ^  s±L,nm 


-I  s 

jeNNf  L',rim' 


Xm'^Z'0^1  +  A'l  +  X  +  AL',n'm’ 


k&NN: 


c 


L',n'm',ij 

L.nm 


(18) 


_  y  y  bj,  ,  ,cL'’n'm'‘i  -  y  y  f  ,  ,f 

/  ;  /  ,  u  L,rimK~'  L,nm  rim!  J  L 


nm  ,i 

nm'J  L.nm 


j£NNj  L',n'm'  n'= 0  m'=-n 

Now  B‘iJnm  is  defined  as  the  solution  to  (8)  when  there  are  only  two  neighboring  particles 

(other  particles  are  not  present)  and  their  displacement  is  purely  tensile  with  no  rotation 
or  far  field.  Thus, 


BiJ  =  -  y  (bj’‘  +  r,’J’iS  8  )cL’n'm’y 

^ L,nm  V  L',n'm'  ^  /  m'  UL'QUn'\  J^L,nm  ’ 


\ 

L',n'm' 


(19) 


where  the  tensile  component  of  displacement  is  given  by 


= 


Jl 


■(rJ  -/) 


2rjrrji  j 


r.. 

p 


(20) 


with  r.  =  r,  -  rt  the  vector  pointing  from  the  center  of  particle  i  to  the  center  of  particle  j. 


We  can  substitute  (19)  into  the  left  hand  side  of  (18)  and  solve  for  A[  nm  as 
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a:  ....  = 


II 

jeNNj  L',n'm' 


{/Jm’  A'l  +  +  X  >W 


keNNj 

k*i 


c 


L',n'm',ij 

L.nm 


(21) 


X  1  X  1  d  j  s~i  L'  ,rim  ,ij  X  '  X  1  j-»  rrirri  ,i 

/  j  /  j  ^ L' .rim  ^ L.nm  /  j  /  j  rim  J  L,nm 


j&NNj  U y rim'  n'= 0  m'=—ri 

Subtracting  out  the  tensile  component  from  the  displacement  in  (21)  allows  the  use  of 
much  lower  values  of  L  because  only  shear  displacements  remain  for  the  nearest 
neighbors.  In  the  computer  code,  we  solve  (19)  to  high  order  in  L.  This  is  facilitated  by 
solving  (19)  for  a  unit  displacement  in  a  reference  frame  where  the  two  particles  are 
aligned  along  the  z  axis  so  that  there  is  cylindrical  symmetry  and  only  m  =  0  values  are 
required  in  the  multipole  expansion,  thus  reducing  the  number  of  multipole  terms  from 
O (Lmax2)  to  O (Lmax).  The  B\'’nm  are  then  derived  from  the  calculated  B‘^n0  by  applying  the 

proper  rotation  using  well-known  rotation  properties  of  the  vector  spherical  hannonics 
(Varshalovich  et  ah,  1988)  and  scaling  according  to  the  actual  displacement  of  the 
particles.  The  coefficients  A\  nm  are  calculated  from  (21)  through  an  iterative  procedure 


and  used  in  (17)  to  calculate  the  multipole  coefficients  B[  nm  . 

The  nearest  neighbor  treatment  allows  the  use  of  a  much  smaller  Lmax  for  an 
accurate  treatment  of  a  particle  pack  and  hence  greatly  reduces  the  memory  requirements 
and  computational  time  for  each  iteration.  However,  the  number  of  iterations  required  to 
reach  convergence  is  increased.  In  a  later  section  we  compare  the  treatment  of  particle 
packs  with  and  without  the  nearest  neighbor  treatment  and  demonstrate  the  computational 
advantages  of  using  this  efficiency  measure  for  highly  filled  packs. 
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5.  Infinite  particle  packs  consisting  of  repeating  cells  of  particles 

Up  to  this  point  we  have  discussed  the  theory  for  a  system  consisting  of  a  finite 
number  of  particles  N.  The  PARMECH  code  also  treats  pseudo-infinite  systems  where 
the  particle  configuration  is  made  by  stacking  a  cell  into  a  three-dimensional  grid  of 
rectangular  unit  cells  consisting  of  N particles.  The  centers  of  all  the  particles  must  fit 
inside  the  unit  cell  in  such  a  way  that  particles  do  not  overlap  with  particles  of  the 
neighboring  cell.  A  slice  through  two  neighboring  unit  cells  of  particles  is  shown  in 
Appendix  B. 

The  formalism  remains  as  described  above  with  the  contribution  from  other  cells 
included  in  the  off-center  expansion  coefficients.  Thus,  the  “pseudo-infinite”  system  is 
actually  a  very  large  finite  system  made  up  of  the  repeated  unit  cells  that  consists  of  an 
inner  region  where  the  relative  displacement  fields  in  each  of  the  cells  is  identical.  This 
inner  region  is  much  larger  than  the  dimension  of  the  unit  cell. 

Only  the  multipole  coefficients  B\  nm  for  the  unit  cell  need  be  treated,  but  the  off- 

center  expansion  coefficients  C.\  Jj  must  include  contributions  from  each  of  the 

neighboring  cells.  Appendix  B  describes  an  efficient  technique  to  calculate  the  off-center 
expansion  coefficients. 

6.  Mechanical  properties  of  a  collection  of  particles 

The  code  will  estimate  the  mechanical  properties  of  a  random  finite  collection  of 
particles  in  an  approximately  spherical  geometry.  Given  a  random  spherical  finite  pack 
of  particles  and  uniaxial  strains  in  the  elastic  medium  far  from  the  particles,  the  code 
calculates  the  displacement  of  each  particle  in  the  pack  as  well  as  estimates  an  effective 
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modulus  and  Poisson’s  ratio  for  the  pack.  This  assumes  the  pack  is  isotropic.  The 
mechanical  properties  of  the  pack  are  detennined  by  a  fit  to  the  exact  solution  of  a 
spherical  body  embedded  in  an  elastic  medium  (Goodier,  1933).  Details  are  given  in 
Appendix  C. 

For  an  infinite  pack,  the  mechanical  properties  of  the  system  are  determined  from 
the  average  stress  and  strain  tensor  seen  in  the  unit  cell,  and  are  derived  in  Appendix  D. 

7.  Method  of  solution  of  the  linear  set  of  equations 

We  employ  an  iterative  method  of  solution  for  the  linear  set  of  equations  (13).  An 
initial  guess  is  provided  for  all  B\  nm ,  that  guess  being  those  multipole  moments 
appropriate  for  an  isolated  sphere  in  an  infinite  medium  with  uniform  applied  fields.  The 
B'l  nm  for  each  particle  i  is  then  recalculated  by  setting  its  value  equal  to  the  right  hand 

side  of  (13).  The  convergence  can  be  surprisingly  fast  depending  on  how  close  the 
particles  are.  For  example,  a  13000-particle  pack  with  a  volume  concentration  of 
particles  C  =  0. 10  will  converge  in  16  iterations  while  a  large  pack  with  volume 
concentration  of  C  =  0.50  may  take  more  than  100  iterations  to  converge. 

8.  Results 

In  this  section  we  summarize  some  of  the  simulations  that  have  been  done  with 
the  PARMECH  code.  We  begin  by  comparing  the  code  results  to  exact  calculations  of 
model  systems  that  have  been  studied  in  the  literature.  This  includes  forces  on  two 
spherical  particles  and  systems  made  up  of  spherical  particles  arranged  as  infinite  cubic 
lattices. 
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We  also  present  results  for  random  arrangements  of  monomodal  and  bimodal 


packs  of  particles.  Care  was  taken  in  the  construction  of  dense  packs  to  prevent  building 
into  the  packs  order  that  would  result  in  artificial  crystallinity.  The  packs  were  generated 
with  a  method  to  build  a  pack  of  touching  particles  (Webb  et  al.,  2006).  The  pack  was 
then  allowed  to  expand  into  a  fixed  volume,  keeping  all  particle  sizes  constant,  using  a 
random  ballistics/collision  dynamic  simulation.  The  dynamic  randomization  was  run 
long  enough  to  allow  the  total  distance  of  travel  of  a  particle  to  be  many  times  the 
distance  across  the  confining  volume.  The  result  is  a  detailed  knowledge  of  the  material 
microgeometry,  with  a  user-specified  volume  fraction  (always  less  than  the  packing 
fraction),  consistent  with  fully  random  fabrication  of  the  material. 

8. 1  Simple  two-particle  system 

To  study  the  convergence  of  the  multipole  solutions  and  the  accuracy  obtained  for 
a  given  Lmax,  we  compared  our  solutions  with  exact  two-sphere  solutions  derived  in 
bispherical  coordinates  by  Shelley  and  Yu  (1966).  Table  I  shows  the  excellent  agreement 
of  the  PARMECH  predictions  with  the  exact  bispherical  solutions.  The  parameter  s, 
mentioned  earlier,  is  a  dimensionless  variable  that  indicates  the  proximity  of  the  spheres 
and  is  defined  as 

e  =  (surface-to-surface  distance)/(diameter  of  smallest  sphere). 

As  expected,  when  the  spheres  are  close  together,  much  higher  Lmax  is  required  to  achieve 
high  accuracy,  hence  the  importance  of  the  nearest  neighbors  treatment  discussed  in 
previous  sections. 


(22) 
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Table  I.  Comparison  of  PARMECH  Predictions  with  Exact  Results  for  the  Forces 
between  Two  Equal  Sized  Spheres  With  Equal  and  Opposite  Uniaxial 
Displacements 


Matrix  Poisson  Ratio  a  =  0.5 


Particle  Proximity  e 

Exact  Force 

PARMECH 
Calculated  Force 

Error 

(%) 

Lmax 

1.0 

3009.6747 

3009.6747 

0.0 

4 

0.1 

8718.9434 

8686.6202 

0.4 

8 

8718.8904 

0.0006 

16 

8718.9417 

0.00002 

26 

0.01 

52993.6830 

47477. 

10. 

21 

50694.3 

4. 

26 

52982.554 

0.02 

50 

52993.6830 

0.0 

100 

0.001 

479048.65 

474229.41 

1.0 

200 

479048.55 

0.00002 

400 

8.2  Infinite  cubic  lattice  arrangements  of  particles 

In  order  to  provide  a  rigorous  testing  of  the  code,  we  have  compared  the  code 
results  to  calculations  in  the  literature  for  infinite  lattice  arrangements  of  SC,  BCC,  and 
FCC  symmetry  (Nunan  et  al.,  1984;  Sangani  et  ah,  1987).  For  these  cubic  lattices  the 
effective  elasticity  tensor  can  be  written  as 

Cijki  =  (A+juf)SjjSki  +  /u(l  +/3){5ik5ji  +  diiSjk)  +  ju(a-  ffifki , 
where  A  and  //  are  the  Lame  constants  of  the  elastic  medium,  Syu  is  unity  if  all  the 
subscripts  are  equal  and  zero  otherwise,  and  a,  J3,  and  /depend  on  the  lattice  arrangement 
and  elastic  medium  Poisson  ratio  a. 

The  PARMECH  code  was  used  to  calculate  a,  fl  and  /for  different  volume 
fractions  C  and  Poisson’s  ratio  a  for  the  simple  cubic  (SC),  body-centered  cubic  (BCC), 
and  face-centered  cubic  (FCC)  infinite  lattices.  A  comparison  of  code  results  to  literature 
values  is  given  in  Tables  II-IV.  The  literature  results  are  read  off  of  small  graphs  and  are 
thus  not  accurate  at  the  third  decimal  place. 


(23) 
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The  excellent  agreement  of  PARMECH  with  literature  values  provides 


confidence  that  the  PARMECH  code  will  predict  correct  results  for  random  packs.  In 
generating  the  effective  elasticity  tensor  for  these  tables,  a  complex  far  field  strain  (i.e., 
Exx  =  0.03,  Eyy  =  -0.06,  Ezz  =  0.1,  Exy  =  0.05,  Exz  =  0.1,  Eyz  =  0.2)  was  applied  to  the 
pack  in  the  PARMECH  code  simulation  so  that  non-trivial  shear  and  tensile  stresses  were 
applied  between  particles,  thus  ensuring  that  the  correct  solutions  was  not  the  result  of  the 
problem  symmetry. 


Table  II.  Comparison  of  PARMECH  Results  with  Results  of  Nunan  et  al.,  (1984) 


and  Sangani  et  al.,  (1987)  for  the  Simple  Cubic  (SC)  Infinite  Lattice 


Volume  fraction  C 
Poisson  ratio  o 

a 

P 

1 

C=0.4 

PARMECH 

2.69 

1.25 

a=0.3 

NK 

2.71 

1.27 

C=0.4 

PARMECH 

4.43 

1.29 

4.73 

a=0.45 

NK 

4.35 

1.27 

4.76 

C=0. 47787 

PARMECH 

5.12 

2.06 

Q 

ll 

o 

UJ 

NK 

5.25 

2.01 

C=0. 47787 

PARMECH 

9.98 

2.10 

5.62 

a=0.45 

NK 

10.4 

2.01 

5.67 

SL 

10.0 

2.10 

5.53 

Table  III.  Comparison  of  PARMECH  Results  with  Results  of  Nunan  et  al.,  (1984) 

|  for  the  Face  Centered  Cubic  (FCC) 

I  Infinite  Lattice 

Volume  fraction  C 
Poisson  ratio  a 

a 

P 

y 

C=0.4 

PARMECH 

1.32 

1.66 

1.46 

a=0.3 

NK 

1.39 

na* 

1.48 

C=0.4 

PARMECH 

1.50 

2.13 

6.35 

a=0.45 

NK 

1.53 

na* 

6.33 

C=0.6 

PARMECH 

4.61 

3.45 

o=0.3 

NK 

3.48 

na* 

3.37 

C=0.6 

PARMECH 

4.48 

7.30 

14.2 

o=0.45 

NK 

4.53 

na* 

na*  Nunan-Keller  results  for  [5  are  incorrect  according  to  Sangani  (1987). 
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Table  IV.  Comparison  of  PARMECH  Results  with  Results  of  Nunan  et  al.,  (1984) 
for  the  Body  Centered  Cubic  (BCC)  Infinite  Lattice _ _ 


Volume  fraction  C 
Poisson  ratio  a 

a 

P 

y 

C=0.4 

PARMECH 

1.30 

1.69 

1.49 

a=0.3 

NK 

1.32 

1.62 

1.45 

C=0.4 

PARMECH 

1.46 

2.18 

6.39 

a=0.45 

NK 

1.48 

2.11 

6.40 

C=0.55 

PARMECH 

2.55 

3.76 

2.89 

a=0.3 

NK 

2.53 

3.71 

2.86 

C=0.55 

PARMECH 

2.97 

5.77 

a=0.45 

NK 

2.98 

5.41 

8.3  Random  monomodal  particle  packs 

PARMECH  simulations  have  been  performed  on  random  particle  packs  with  both 
a  finite  spherical  geometry  and  infinite  packs  made  up  of  repeating  unit  cells.  Results 
for  monomodal  packs  are  summarized  in  Table  V.  The  Poisson  ratio  in  these  calculations 
is  0.5.  We  list  the  order  Lmax  used  in  the  general  multipole  expansion  and  the  multipole 
expansion  for  the  nearest-neighbors  treatment.  The  average  minimum  s  (see  (22))  listed 
in  this  table  is  the  average  over  all  particles  in  the  pack  of  the  minimum  s  for  each 
particle  with  its  closest  neighbor. 

The  PARMECH  simulations  of  the  finite  packs  in  Table  V  used  a  uniaxial  tensile 
strain  as  described  in  Appendix  C.  The  simulations  of  the  infinite  pack  applied  a 
complex  far-field  strain  (i.e.,  Exx  =  0.03,  Eyy  =  -0.13,  Ezz  =  0.1,  Exy  =  0.05,  Exz  =  0.1, 
Eyz  =  0.2)  to  the  pack.  The  average  stresses  and  strains  in  the  medium  were  calculated  as 
described  in  Appendix  D  for  the  infinite  medium.  For  an  isotropic  homogeneous  linear 
elastic  medium,  the  average  stress  is 
related  to  the  average  strain  by 
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^=T7^(Sf«>5»  +  2^- 

1  +  aeff  k 

The  simulations  were  performed  with  the  condition  ^  fkk  =  0 ,  so  that 

k 

fr 

Tj  =  2 neffE“j  or  /ueff  =  •  Because  the  packs  do  not  represent  a  perfectly  isotropic 

medium,  we  can  obtain  an  independent  effective  shear  modulus  for  each  of  the  six 
independent  values  of  ij.  The  values  shown  in  Table  V  are  an  average  of  the  six  values 
with  their  standard  deviations. 

For  the  pack  with  C  =  0.06  and  C  =  0. 1,  a  nearest  neighbors  treatment  is  not 
needed  and  Lmax  =  8  is  adequate  for  good  accuracy  because  the  average  minimum  s  for 
these  packs  is  greater  than  0. 1,  a  particle  separation  that  was  treated  with  good  accuracy 
with  Lmax  =  8  in  Table  I.  A  comparison  of  results  with  identical  packs  for  C  =  0.3  and  C 
=  0.4  (packs  with  71  and  49  particles,  respectively)  indicates  that  a  maximum  value  of  L 
=  8  for  these  cases  is  probably  adequate  to  get  the  needed  accuracy  when  the  nearest- 
neighbors  method  is  used.  However,  improved  accuracy  is  achieved  for  C  =  0.5  with  L  = 
12,  as  seen  for  the  pack  with  50  particles. 

The  degree  that  the  pack  deviates  from  isotropy  can  be  estimated  by  the  standard 
deviation  of  the  effective  modulus.  It  can  be  seen  that  the  packs  with  the  higher  volumes 
of  solids  have  the  highest  deviations  from  isotropy.  Also,  as  expected,  the  packs  with 
larger  numbers  of  particles  generally  have  a  lower  standard  deviation. 
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Table  V.  PARMECH  Predictions  of  Effective  Shear  Modulus  Divided  by  Matrix  Shear 
Modulus  (  jueff  /  jUmatrix  )  for  Monomodal  Packs 


Poisson’s  Ratio  cr=  0.5 


Volume  fraction 
C  of  Pack 

Pack 

Geometry 

Average 

Minimum 

8 

L  max 

Nearest 

Neighbors 

h  mux 

Number 

of 

Particles 

Meff  /  M  matrix 

C=  0.06 

Spherical 

0.67 

4 

no 

1279 

1.18 

Spherical 

0.67 

4 

no 

12820 

1.18 

Infinite 

0.53 

8 

200 

30 

1.16±0.01 

Infinite 

0.36 

8 

200 

147 

1.17±0.002 

C=0.1 

Spherical 

0.50 

4 

200 

1278 

1.28 

Spherical 

0.50 

8 

200 

12820 

1.28 

Infinite 

0.26 

8 

200 

39 

1.30±0.02 

Infinite 

0.23 

8 

200 

148 

1.32±0.02 

C=  0.2 

Infinite 

8 

48 

1.80+0.13 

Infinite 

145 

1.75+0.05 

C=  0.3 

Infinite 

0.039 

8 

200 

71 

2.51+0.38 

Infinite 

0.039 

12 

200 

71 

2.51+0.38 

Infinite 

0.052 

10 

200 

161 

2.72+0.19 

C=  0.4 

Infinite 

0.029 

8 

200 

49 

4.05+0.38 

Infinite 

0.029 

12 

200 

49 

4.06+0.37 

Infinite 

0.027 

12 

200 

156 

4.47+0.55 

Infinite 

0.027 

12 

200 

312 

4.34+0.10 

C=  0.5 

Infinite 

0.0083 

8 

200 

50 

8.66+1.15 

Infinite 

0.0083 

12 

200 

50 

8.74+1.10 

Infinite 

0.0121 

12 

200 

167 

7.80+2.15 

Infinite 

0.0109 

12 

200 

318 

7.76+0.54 

C=  0.55 

Infinite 

12 

47 

13.28+1.09 

Infinite 

12 

200 

175 

18.08+5.61 

Infinite 

12 

200 

332 

12.34+1.27 

8.4  Accuracy  of  nearest  neighbors  treatment 

As  discussed  previously,  the  nearest-neighbors  approximation  treats  the  tensile 
displacements  of  particles  in  close  proximity  with  a  high  order  of  multipole  terms,  and  is 
most  useful  when  the  Poisson’s  ratio  is  close  to  Vi.  We  compare  the  results  of  the  code 
using  the  nearest-neighbors  treatment  to  results  obtained  to  high  order  without  the 
nearest-neighbors  approximation.  Table  VI  shows  the  results  of  these  comparisons  with 


26 


monomodal  particle  packs  at  0.50  and  0.55  volume  fraction  of  solids.  Again,  the 
Poisson’s  ratio  is  Vi. 


Table  VI. 

Comparison  of  Nearest  Neighbors  Treatment  to  High  Order  Treatment  of  Highly 

Filled  Monomodal  Packs 

Poisson’s  Ratio  <j=  0.5 

Aver- 

Near- 

Memory 

Volume 

age 

est 

for 

Number  of 

Rela- 

fraction 

Min- 

Number 

Neigh- 

Matrix 

Iterations 

tive 

E/E  matrix 

Cof 

imiim 

of  Part- 

bors 

Storage 

to  Conver- 

CPU 

Pack 

e 

icles 

Lmax 

Lmax 

(Gb) 

gence 

Time 

O 

II 

o 

bi 

0.0083 

50 

12 

200 

0.85 

428 

1.0 

8.742±1.10 

12 

no 

0.83 

89 

0.23 

8.1 1+0.82 

16 

no 

1.87 

121 

0.56 

8.43±0.92 

20 

no 

3.68 

152 

1.35 

8.58±0.97 

24 

no 

6.54 

181 

2.85 

8.67±1.00 

28 

no 

10.83 

209 

4.66 

8.717+1.03 

32 

no 

16.85 

234 

9.13 

8.744±1.04 

36 

no 

24.95 

258 

15.7 

8.761+1.06 

C=  0.55 

0.0047 

47 

12 

200 

0.84 

(in  RAM) 

416 

1.0 

13.29+1.08 

12 

no 

0.82 

(in  RAM) 

97 

0.41 

10.59+1.71 

16 

no 

1.88 

133 

2.26 

1 1.43+1.87 

20 

no 

3.75 

167 

5.13 

1 1.97+1.86 

24 

no 

6.78 

199 

10.9 

12.33+1.78 

28 

no 

11.38 

233 

22.1 

12.57+1.68 

32 

no 

17.96 

266 

40.9 

12.75+1.59 

36 

no 

27.17 

295 

65.3 

12.87+1.52 

From  Table  VI  it  can  be  seen  that  the  nearest-neighbors  treatment  is  effective  at 
predicting  mechanical  properties  of  highly  filled  packs  using  a  relatively  low  order  in  the 
multipole  expansions  with  less  memory  and  CPU  requirements  than  the  simulations 
without  the  nearest-neighbors  option.  It  is  also  apparent  that  Lmax  =  12  is  adequate  to 
achieve  good  accuracy.  We  note  that  the  predicted  results  with  the  nearest  neighbors 
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treatment  is  not  expected  to  be  exactly  the  same  as  a  high  order  multipole  simulation 
without  the  nearest  neighbors  treatment  since  the  nearest  neighbor  treatment  only  uses  a 
high  order  multipole  on  the  dominant  tensile  displacement,  and  does  not  use  a  higher 
order  on  the  shear  displacements  between  particles. 

8.5  Random  bimodal  particle  packs 

One  of  the  reasons  the  PARMECH  code  was  developed  was  to  help  guide  the 
construction  of  constitutive  models  that  might  use  statistical  averages  for  the  arrangement 
of  particles  in  a  material  to  estimate  its  mechanical  properties.  Such  a  task  can  become 
increasingly  challenging  when  different  sizes  of  particles  can  interact  in  a  medium. 
Therefore,  we  have  run  PARMECH  simulations  on  a  number  of  bimodal  packs  with 
particle  volume  fractions  from  0.2  to  0.6.  These  have  been  run  with  the  volume  ratio  of 
coarse  to  fine  particles  at  3  and  9.  The  ratio  of  diameters  of  coarse  to  fine  particles  varies 
from  2.5  to  7.5. 

Results  are  summarized  in  Table  VII.  It  will  be  noted  that  no  results  are  given  for 
some  systems  at  the  0.6  volume  fraction.  This  is  because  the  code  sometimes  does  not 
converge  for  high  volume  fractions  when  the  nearest-neighbors  treatment  is  used.  It 
appears  that  the  stiff  interactions  of  many  particles  in  close  proximity  sometimes  results 
in  particle  displacements  that  grow  during  the  iteration  process.  We  have  been  able  to 
improve  the  convergence  behavior  by  switching  on  the  nearest-neighbors  treatment 
gradually,  but  there  are  still  occasions  when  the  code  does  not  converge. 

We  have  compared  the  calculations  to  the  predictions  of  an  effective  medium 
theory  where  the  large  particles  are  assumed  to  be  in  an  effective  homogeneous  medium 
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with  mechanical  properties  of  the  small  particles  in  the  elastic  matrix.  These 
comparisons  are  shown  in  Table  VII. 

In  the  table  the  volume  fraction  C  of  the  pack  is  given  as  the  sum  of  the  coarse 
and  fine  volume  fractions  of  particles  Cc  and  C/.  The  modulus  of  the  effective  medium  is 
calculated  from  a  simulation  of  a  monomodal  pack  at  a  volume  fraction  C^that  is  given 
by  the  fine  particles  filling  a  volume  of  the  pack  that  consists  of  the  total  volume  minus 
the  volume  of  coarse  particles  Ceff=  C/7 f  1  -  Cc).  The  calculation  of  the  modulus  of  the 
coarse  particles  was  made  using  the  same  particle  pack  as  the  bimodal  system,  but  with 
the  fine  particles  removed.  Clearly,  when  the  particles  sizes  differ  by  a  large  amount, 
such  an  effective  medium  approximation  should  be  accurate.  However,  when  the 
surface-to-surface  distances  of  the  larger  particles  are  on  the  order  of  the  size  of  the 
smaller  particles,  then  the  critical  high-stress  region  between  particles  cannot  be  treated 
as  an  effective  medium.  For  this  case,  the  effective  medium  theory  should  always  give  a 
larger  estimate  for  the  modulus  since  the  region  between  particles  is  assumed  to  have  the 
larger  modulus  of  the  effective  medium  of  the  fine  particles.  That  is  what  is  seen  in  these 
calculations.  When  one  compares  the  average  minimum  surface-to-surface  distances 
represented  by  the  non-dimensional  parameter  s  to  the  particle  size  of  the  small  particles, 
it  is  somewhat  surprising  that  the  effective  medium  theory  does  so  well.  Of  coarse,  the 
agreement  is  worse  as  the  volume  fraction  of  particles  increases. 

It  is  our  plan  to  use  these  and  additional  calculations  of  the  PARMECH  code  to 
assist  in  developing  a  constitutive  theory  based  on  statistical  pack  properties  for  the 
effective  modulus  of  composite  materials  containing  rigid  spherical  particles  of  different 
sizes. 
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Table  VII.  PARMECH  (PM)  Calculations  of  Modulus  for  Bimodal  Particle  Packs  with  Matrix  Poisson’s  Ratio  Equal  to  1/2 

Comparisons  to  an  Effective  Medium  Theory  (EFT) 

Bimodal  Pack 

Coarse  Particles  Alone 

(Pack  with  fines  removed) 

Effective  medium, 
202  fine  particles 

C/Cf 

C 

Rc/Rf 

Number 
of  Coarse 
Particles 

Number 
of  Fine 
Particles 

PM 

E/E0 

EFT 

E/EeffxEef/E0 

Percent 

Differ¬ 

ence 

Cc 

Average 

Minimum 

£ 

PM 

E/Eeff 

Ceff 

PM 

Eef/Eo 

3 

0.2 

2.5 

54 

281 

1.897 

1.916 

1.0 

0.15 

0.117 

1.636 

0.0588 

1.166 

3 

0.2 

5.0 

19 

791 

1.755 

1.781 

1.5 

0.15 

0.130 

1.521 

0.0588 

1.166 

3 

0.2 

7.5 

10 

1406 

1.748 

1.767 

1.1 

0.15 

0.173 

1.509 

0.0588 

1.166 

3 

0.4 

2.5 

54 

281 

4.700 

4.562 

-2.9 

0.3 

0.022 

3.064 

0.1429 

1.499 

3 

0.4 

5.0 

19 

791 

4.074 

4.203 

3.2 

0.3 

0.045 

2.823 

0.1429 

1.499 

3 

0.4 

7.5 

10 

1406 

4.547 

5.127 

12.7 

0.3 

0.019 

3.443 

0.1429 

1.499 

3 

0.5 

2.5 

54 

281 

8.341 

8.670 

3.9 

0.375 

0.015 

4.904 

0.2 

1.766 

3 

0.5 

5.0 

19 

791 

7.982 

9.026 

13.1 

0.375 

0.006 

5.105 

0.2 

1.766 

3 

0.5 

7.5 

10 

1406 

9.489 

11.122 

17.2 

0.375 

0.010 

6.291 

0.2 

1.766 

3 

0.6 

2.5 

54 

281 

nc 

15.693 

nc 

0.45 

0.006 

7.012 

0.2727 

2.279 

3 

0.6 

5.0 

19 

791 

21.543 

24.909 

15.6 

0.45 

0.005 

11.13 

0.2727 

2.279 

3 

0.6 

7.5 

10 

1406 

nc 

17.539 

nc 

0.45 

0.002 

7.837 

0.2727 

2.279 

9 

0.2 

2.5 

68 

118 

1.888 

1.874 

-0.8 

0.18 

0.106 

1.761 

0.0244 

1.065 

9 

0.2 

5.0 

39 

542 

1.761 

1.770 

0.5 

0.18 

0.134 

1.664 

0.0244 

1.065 

9 

0.2 

7.5 

15 

703 

1.774 

1.778 

0.2 

0.18 

0.096 

1.671 

0.0244 

1.065 

9 

0.4 

2.5 

68 

118 

4.133 

4.203 

1.7 

0.36 

0.024 

3.571 

0.0625 

1.176 

9 

0.4 

5.0 

39 

542 

4.471 

4.588 

2.6 

0.36 

0.019 

3.898 

0.0625 

1.176 

9 

0.4 

7.5 

15 

703 

3.834 

4.015 

4.7 

0.36 

0.028 

3.411 

0.0625 

1.176 

9 

0.5 

2.5 

68 

118 

9.261 

9.742 

5.2 

0.45 

0.011 

7.629 

0.0909 

1.269 

9 

0.5 

5.0 

39 

542 

6.480 

6.999 

8.0 

0.45 

0.014 

5.481 

0.0909 

1.269 

9 

0.5 

7.5 

15 

703 

7.598 

8.318 

9.5 

0.45 

0.010 

6.514 

0.0909 

1.269 

9 

0.6 

2.5 

68 

118 

nc 

19.742 

nc 

0.54 

0.005 

13.71 

0.1304 

1.445 

9 

0.6 

5.0 

39 

542 

18.363 

22.277 

21.3 

0.54 

0.004 

15.47 

0.1304 

1.445 

9 

0.6 

7.5 

15 

703 

nc 

21.701 

nc 

0.54 

0.001 

15.07 

0.1304 

1.445 

nc  =  no  convergence 
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9.  Summary 

The  PARMECH  computer  code  calculates  the  stress  and  strain  fields  of  a 
deformed  medium  consisting  of  a  pack  of  rigid  spherical  particles  embedded  in  an  elastic 
matrix.  The  desired  stress  and  strain  can  be  calculated  at  any  point  in  the  medium  to 
within  a  user-specified  accuracy,  limited  only  by  the  computation  time.  Average 
mechanical  properties  of  the  medium  are  also  output  by  the  code.  The  code  can  be  used 
to  simulate  systems  consisting  of  thousands  of  particles  in  a  finite  pack  or  an  infinite  pack 
consisting  of  a  repeating  rectangular  cell  of  particles. 

A  multipole  expansion  technique  is  used  to  solve  the  equations  of  small 
deformation  for  the  elastic  medium  and  consists  of  truncated  sums  of  complete 
orthogonal  vector  spherical  hannonics  that  are  exact  analytical  solutions  for  a  single 
spherical  particle.  A  solution  is  obtained  by  solving  for  the  coefficients  in  the  multipole 
expansion  from  a  large  set  of  simultaneous  linear  equations.  Techniques  are  used  to 
improve  the  convergence  of  the  solution  when  the  particles  are  in  close  proximity  for 
highly  filled  mediums. 

The  PARMECH  code  results  agree  with  exact  solutions  of  configurations 
consisting  of  a  few  particles  as  well  as  infinite  packs  of  particles  in  body-centered  cubic, 
face-centered  cubic,  and  simple  cubic  lattice  arrangements. 

The  code  has  been  used  to  estimate  mechanical  properties  of  a  variety  of 
monomodal  and  bimodal  particle  packs  of  different  packing  densities.  Eventually,  the 
code  will  be  used  to  help  derive  a  constitutive  theory  for  highly  filled  particulate 
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materials  with  broad  size  distributions  of  particles  that  debond  from  the  matrix  during 
loading. 


Appendix  A.  Translation  coefficients  of  Navier  multipole  solutions  valid  on  the 
surface  of  Sphere  i 

In  this  appendix  we  use  addition  theorems  to  derive  off-center  expansion  coefficients 
of  (11).  The  off-center  expansion  coefficient  is  defined  as 


m . ,ij 
^  Lj  ynimi 


2n  1 

\d(p]d (cos (#, ,<Pi )  •  uLj  njmj  (r  -  j Tj ) 


-1 


(Al) 


where  uL  n  m  (r  -  r. )  is  the  single  particle  solution  for  sphere  j  defined  in  (3b) 


evaluated  on  the  surface  of  sphere  i.  The  global  vector  positions  of  the  two  spheres  are 
denoted  by  r  and  r  .  The  position  of  Sphere  j  relative  to  Sphere  i  is  given  by  Rjt  =  r.  -  r. 
We  also  define  a  vector  emanating  from  Sphere  j  to  a  point  on  the  surface  of  Sphere  i  by 


The  off-center  expansion  will  consist  of  two  parts:  the  harmonic  term  V[xnm{rjj) 


and  the  nonharmonic  term  dn  (ri  -a~)V[V  ■  V‘'x n  m  (f/;)]  for  solutions  external  to  the 
sphere.  The  off-center  multipole  expansion  for  the  harmonic  term  is 

i 

oo  «;+ 1  nt 

L  T.D&Z&J  (A2) 
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\rj‘  J 


=0L,=n,-\  m:=- n 


where  the  expansion  coefficients  (R^  )  are  functions  only  of  the  relative  position 


vector  R  „ .  We  will  derive  them  below. 
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The  off-center  multipole  expansion  for  the  nonharmonic  term  is 


(211, -!)>>,■  +D  -  -  Y£(0J,,pJl)  = 


oo  «,+l 

Z  Z  EaS'”' A> 


«,=OL=n,-l 


where  we  will  only  calculate  the  off-center  term  for  Lj  =  rij—  1  since  the  nonhannonic 
tenns  are  zero  for  the  other  two  values  of  Lj.  Thus  the  off-center  expansion  coefficient  is 

(R,)-  SLjnrldnEn^  (R,)  (A4) 

where  dn  is  the  same  as  the  dn  in  (3d). 


The  harmonic  term  coefficients  are  given  by  the  integral 


d11;ZZ'(rji)=  jdf>,jd(cose,)  -l  <*,*,)•?£<*»>«> 


and  the  nonharmonic  term  coefficients  by  the  integral 


jjnj  -\;rij  ,mj 


L,n„mi  (Rji)  =  (2nj-l^nj(nJ+l) 


j"  d(pt  j*  d  (cos  ^ 


Expressions  for  scalar  hannonic  function  off-center  expansions,  which  are  needed  in  (A5) 
and  (A6)  for  the  functions  of  (r/(,  A/(, ^/(),  are  given  by  Varshalovich  (1988): 


l47T(2L  +  1)  yyi  .  jy 

1  (2  L)\  h 


|(2Z  +  2/l)!  afaj+l 
(2 A  + 1)!  R^m 


zP  L+A  L  w>  • 

m~R  ~m) 
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Performing  the  integral  of  (A5)  for  the  harmonic  term  of  the  Navier  off-center  expansion, 
one  obtains 


nLi’nJ’mi(R  \  =  r]LJ’nJ’mj 


L,  L.+ 1 

a,  ‘a, 

1  J  .y 
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L,  +L :  ,m :  —m 


(A8) 


where  the  constant  coefficient  is  given  by 
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(A9) 


where  the 


A  A  A 


vm,  m2  m3y 


are  the  Wigner  3j  symbols  which  are  also  defined  in 


Varshalovich  (1988). 

The  nonharmonic  tenn’s  off-center  expansion  is  more  complicated  but  it  only 
exists  when  Lj  =  nj—  1 .  We  can  write  the  nonharmonic  tenn  as  follows: 
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(A10) 


We  have  used  the  law  of  cosines,  rH  =  a]  +  R~n  -  2 aiRji  cos  a  ,  where 
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(All) 


cos  a  =  cos  0 n  cos  6l  +  sin  0  /;.  sin  0j  cos (<pi  -  O  .. ) 

=  cos  0  cos  0  +  —  sin  0  ..  sin  6i  1  +  )  . 

Jl  1  ^  Jl  1  \  ) 


Then  (A6)  becomes 
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(A  12) 


aj  I  /v  (7  2  ci  R  (/  R  ,(I) 

- ^ - - - jA.  cos  0  ..  cos  Q, - sin  0  /Ye  ■"  sin 


a/ 


a 


a; 


a,J?. 


'^-sin0y/e'  "  sin  6;e 


The  first  term  in  the  final  square  brackets  of  (A  12)  has  no  angular  dependence  and  can  be 
evaluated  as  before,  i. e. ,  as  a  harmonic  term.  The  second  term  has  a  cos#  factor.  The 
third  term  has  a  sin  0ieup‘  and  the  fourth  term,  a  sin  6te  Uf>i  .  The  dot  product  of  the  vector 

spherical  harmonics  turns  the  expression  into  a  sum  of  scalar  spherical  harmonic 
products.  Then  the  second,  third,  and  fourth  terms  of  the  final  square  brackets  can  be 
turned  into  pure  spherical  harmonics  via  the  following  recursion  relations  found  again  in 
Varshalovich  (1988). 


cos  0Yhn{0,(p)  =  . 


(/  -  m  + 1)(/  +  m  + 1) 


(2/  +  l)(2/  +  3) 


W  0,<P) 
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(A13a) 
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~i(P  V  (  P  sn\  —  /(^  ^+1)(7  in +  2) 
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(A13c) 


In  using  these  recursion  relations,  we  note  that  in  the  integral  of  (A  12)  the  vector 
spherical  harmonic  we  are  integrating  is  complex  conjugated  and  so  we  use 


(0, ,  <pt )  cos  0,  =  [7^  (0. ,  <pt )  cos  0,  ]* 

(0: ,  <p, )  sin  =  [?£  (0, ,  cpi )  sin  0^*  ]* 
S  (0, ,  P, )  sin  0., fi*  =  [fi  (0, ,  p, )  sin  0,e**  ]* 


(A  14) 


The  entire  integral  of  (A  12)  is  laborious  but  straightforward: 
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(A  15) 


Detailed  evaluation  of  g"Lj’™Jm  shows  that  this  coefficient  is  identically  zero.  It  will  not  be 


written  out  explicitly  in  the  following  equations. 


The  coefficient  e"”mm  of  the  second  line  of  (A  15)  is  defined  by 
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(: 2Li+2n,.)\ 
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(A  16) 


(A  16a) 


(A  16b) 
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e3(n  j,mj;L,nj,mi)  = 
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There  are  two  coefficients  in  the  third  line  of  (A  15).  The  coefficient  /YY  is  given  by 
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The  second  coefficient,  ,  in  the  third  line  of  (A  15)  is  given  by 
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Appendix  B.  Treatment  of  a  pseudo-infinite  particle  pack  via  repeating  cells 

A  pseudo-infinite  particle  pack  is  created  by  making  a  collection  of  particles  with 
centers  that  lie  in  a  cubic  cell  with  volume  Z  and  then  repeating  the  cell  a  large  number 
of  times  in  all  directions.  It  is  assumed  that  Nc  (with  Nc  »1)  cells  from  a  relatively-small 
central  portion  of  this  pack  undergo  identical  distortions  when  the  whole  pack  is 
distorted.  We  refer  to  this  as  a  pseudo-infinite  pack.  This  geometry  allows  us  to  use  the 
same  mathematical  framework  as  with  a  finite  pack  that  has  asymptotic  deformations  far 
from  the  particles. 

The  particle  pack  must  be  made  so  that  a  particle  in  a  cell  does  not  overlap  a 
particle  in  a  neighboring  cell.  Figure  B-l  shows  two  adjacent  cells  from  a  two 
dimensional  slice  through  a  pack. 
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Figure  B-l.  Cross  section  of  two  identical  adjacent  cells  of  particles 


Great  computations  efficiencies  can  be  realized  because  only  the  multipole  coefficients 
for  a  single  cell  need  be  stored  as  each  cell  is  identical  and  because  analytical  expressions 
can  be  used  to  sum  the  effect  of  exterior  cells  on  a  central  cell. 

We  delineate  the  cell  by  the  subscript  s  with  5  =  0  being  the  center  cell,  4  being 
the  ith  particle  in  cell  number  s.  Since  each  cell  is  identical,  the  multipole  coefficients  do 
not  depend  on  the  cell  in  which  they  are  contained  and  B'f  nm  =  B'J  nm  .  With  this 
nomenclature  we  can  rewrite  (11)  for  the  particles  in  the  center  cell  as 
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2  ri 

D i  _  i  c  e  ,  z  o*  c  _  \  ’  \  ’  77  rn'm',i  _  \  1  \  1  n/  s-^L' ,rim' ,iQis 

D L,nm  ~  7 m°LO°n\  Y  m°L\°n\  Zu  Lu  F n'm'  J L ,nm  Zu  /  4  L'.n'm'^ L.nm 

n'= 0  m'=—n'  s^O  L',n'm' 

-zf  =  KA,A,  +  VJL A,  -£ 

5  j=\,j^i  L',n'm'  n'= 0  m'=—n' 


(Bl) 


En*  J^iL',n'm',ii  \  1  \  1  n)  fiL',rim',i 

^ L',n'm,{^L,nm  ~  Zj  Z  >  L'  .rim' ^  L.nm 

j=\Jld  L' ,n'm' 


L'.n'm' 


where 


s~i  L' ,n'm' ,ii  _  \  1  s~i  L'  ,rim'  ,i0is 

^  L,nm  /  v  ^  L,nm 


s*  0 


(B2a) 


and 


C 


L',n'm',ij 

L.nm 


Zsi L' ,n'm  ,i0js 
^ L,nm 
s 


(B2b) 


Equation  (Bl)  is  the  same  form  as  (1 1)  with  the  trivial  difference  that  the  right  hand  side 
of  (B 1)  depends  on  B[  nm  so  that  the  same  solution  technique  can  be  applied  to  both  the 
finite  particle  pack  and  the  infinite  particle  pack  of  repeating  cells. 

To  calculate  the  Qj"’’  one  could  do  a  direct  summation  over  cells  according  to 

(B 1).  This,  however,  can  be  computationally  intensive  for  low  values  of  L  where  a  large 
number  of  cells  contribute.  We  have  chosen  another  method.  This  technique  will  sum 
the  effect  of  all  particles  in  a  cell  (or  some  smaller  portion  of  a  cell  that  we  will  call  a 
subcell)  and,  using  off  center  expansions,  will  write  the  effect  in  tenns  of  a  multipole 
expansion  centered  in  that  cell  (or  subcell).  The  summation  over  all  the  cells  can  then  be 
done  a  single  time,  saving  the  computational  expense  of  summing  the  effect  of  each 
particle  over  all  the  cells. 

We  again  label  each  cell  by  the  subscript  s.  As  will  be  shown,  computational 
efficiencies  are  achieved  by  splitting  up  each  cell  into  subcells  that  we  enumerate  with 
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the  subscript  t  as  shown  in  the  Figure  B-2.  We  align  the  boundaries  of  the  cells  so  that 
particle  i  lies  in  the  center  of  the  center  cell. 


Ce 

l\s 

Particle 

insubc 

M,, 

i\ 

Express  the  effect 
of  particles  in 
subcell  t  in  terms 
of  offcenter 
expansions  at  the 
subcell  center 

7  \ 

Particle  i  at  center 
of  celled 

Figure  B-2.  The  unit  cell  is  divided  into  2x2x2  =  8  subcells  in  the  above  example. 


We  begin  the  derivation  by  rewriting  (9)  for  the  general  solution  to  the  Navier 
equation  with  a  grouping  of  particles  in  their  cells  and  subcells  as 

s(o = £  + Z"j'<2  -*,) 

n=0  m=-n  j  .  , 

In  N  Nsub  _  ' 

= E  fL (^)+2«*(?-**)+Z2  S#*' - K ) 

n=Om=-n  7=1  s* 0  t= 1  yesubcell  t 

where  Nsub  is  the  number  of  subcells  in  a  cell,  N,  is  the  number  of  particles  contained  in 
subcell  t,  and  jSJ  represents  the yth  particle  in  cell  s  and  subcell  1.  The  second  term  on  the 
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right  hand  side  of  (B3)  is  simply  the  sum  of  contributions  from  particles  in  the  center  cell 
and  is  treated  in  the  same  way  as  a  finite  pack.  The  last  tenn  in  (B3)  gives  the 
contribution  from  the  particles  of  all  other  cells.  Its  treatment  is  described  below. 


We  now  write  the  contribution  of  all  the  particles  in  a  subcell  in  terms  of  their  off- 
center  expansion  around  the  center  of  the  subcell,  valid  for  a  point  outside  of  the  subcell. 
The  appropriate  off-center  expansion  is  given  in  Appendix  E.  Using  (3)  and  the  off- 
center  expansions  we  write  the  last  tenn  of  (B3)  as 

N„lb  _ 

XX  X"'"'<r  tf/J-XX  X  YjBkn'm’U,’,n'Jr  -  RjJ 

5^0  t= 1  yesubcell  t  s* 0  t=\  yesubcell  t  L',n',m' 


ZZ  Z  ZBw 

5^0  t= 1  yesubcell  t  L',n',m' 


(B4) 


Z  few  (A  prJ~  WAS  (*„ ,  .<•„  )te 


where  RstJ  is  the  vector  from  the  center  of  the  subcell  labeled  by  s,t  to  the  center  of 
particle j  and  fst  =[rst,Ost,(pst)  is  a  vector  in  spherical  polar  coordinates  centered  at  the 
center  of  the  cell  labeled  by  s,t. 

An  inspection  of  D\\%(Rs  tprs  t)  and  in  Appendix  E  shows 

that  terms  in  (B4)  are  proportional  to  either  -X-  or  4XX'(0sV,fev)  . 


Now  we  use  off-center  expansions  derived  in  Appendix  A  to  expand  these  terms 
at  the  surface  of  a  particle  i  located  at  the  center  of  the  center  cell  (see  Figure  B-2). 


L"+ 1 


=  y  D^"'n"m"YL  (0.,<p.) 

L,nm  nm  \  i  ‘  /  /  / 


(B5) 


Surface  of  i 


L,n,m 
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where  Rs  ti  is  the  vector  from  the  center  of  particle  i  to  the  center  of  the  subcell  labeled 


by  s,t.  The  coefficients  d 


'  L"  ,n"m" 
L,nm  ‘ 


,  and  are  given  in  Appendix  A. 


Substitution  of  (B5),  (B6),  (B7),  and  (B8)  into  (B4),  and  (B4)  into  (B3),  yields  an 
equation  that  describes  the  displacement  field  in  terms  of  vector  spherical  hannonics  in 
the  reference  frame  of  a  particle  i  located  at  the  center  of  the  center  cell.  Applying  the 
boundary  condition  (10)  about  particle  i  and  the  orthonormality  of  the  vector  spherical 
hannonics  (6)  yields  an  equation  with  the  form  of  (Bl)  with 
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5^0  t= 1  L" ,n" ,m" 


(B9) 


-<Av-,<A7'*‘  8SvW-.-(e.«.®«)tev(i.) 

+ 1 1 1  [-  rfA«  <2»'  -  iw^F^r<r%sw-^  (®„/.  ®,, , ) 

5^0  £=1  «",m" 

_  j  o  n' pn"-n' rfn'm'  y  uhI  (T)  I  Z7  /  r>  \ 

^  n'^L’n'-l^  j  ^ s,t  j  J n"-l,n"m"*  n"-n',m'-mn  \  ^ s,t  j 9  s,t  j  /J^ Z,,«m  / 

.  (1  —  A  X^L'>n'm,jj  I  nL',rim',ij 
"r  V1  UijryL,nm  ^  HL,nm 

where  the  second  to  last  term  in  (B9)  is  the  off-center  expansion  contribution  of  sphere  j 
on  sphere  i  in  the  center  cell  and  ft'  nnmm  ,ij  is  a  boundary  correction  due  to  the  fact  that 

sphere  i  is  not  actually  at  the  center  of  the  center  cell. 

It  is  now  evident  why  we  split  a  cell  into  subcells  for  the  treatment  of  off-center 


R 


L" 
s,t  j 


expansions.  Terms  in  (B9)  are  of  the  form  SL„J  and  become  smaller  (much  less  than 

R, , , 


one)  as  the  number  of  subcells  increase.  Therefore  the  sum  can  be  tenninated  at  lower 
values  of  L"  when  a  cell  is  broken  into  a  number  of  subcells. 

Since  (B9)  depends  only  on  the  geometry  of  the  pack,  it  can  be  evaluated  a  single 
time  at  the  beginning  of  a  simulation.  An  iterative  process  is  used  to  solve  for  the 
multipole  coefficients  in  (Bl). 

In  the  method  described  above,  the  contribution  of  all  particles  from  the  center 
cell  are  summed  separately  from  the  contribution  of  all  cells  other  than  the  center  cell. 

We  have  also  found  it  efficient  to  sum  all  the  particles  explicitly  from  the  27  cells  nearest 
the  center,  and  then  sum  the  contribution  of  all  other  cells.  In  that  case,  there  is  no  need 
to  divide  the  cell  into  subcells. 
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Correction  due  to  assumption  that  each  particle  is  at  the  center  of  the  pack 

At  first  glance  it  would  appear  that  there  could  be  convergence  problems  of  the 

sum  over  cells  of  (B9)  due  to  the  fact  that  some  terms  in  the  summation  for  C'ff'  and 


Q,"' ,ij  are  proportional  to 


and  the  number  of  summations  at  a  given  distance 


increases  as  R2s  ti .  This  is  not  the  case  because  the  angular  dependence  of  these  terms 

results  in  cancellation  of  contributions  from  the  symmetry  of  the  cells  around  the  central 
particle.  However,  each  particle  is  not  actually  at  the  central  point  of  the  central 
collection  of  N ]  cubes.  Therefore  there  is  an  additional  contribution  ( (fff  'J  in  (B9)) 


for  terms  that  are  proportional  to 


because  there  is  not  a  cancellation  due  to  the 


symmetry  of  the  pack.  This  extra  term  is  derived  below. 

Let  Rf  be  the  vector  from  the  center  of  the  central  cell  to  particle  i.  The  center  of 

the  pack  considered  in  (B9)  is  therefore  shifted  by  Rt  to  the  true  center  of  the  pack. 

When  evaluating  (B9)  we  are  summing  contributions  from  additional  particles  on  the 
boundary.  For  example,  if  the  z  component  of  Rt  is  a  positive  Rz . ,  then  the  boundary  at  a 
large  +z  direction  contains  a  group  of  particles  defined  in  the  central  cell  by 
S~  =  {particles  j,  such  that  Rz .  -  R,  .  >  L/2  }  . 

The  contribution  from  these  particles  must  be  subtracted  from  (B9).  Similarly,  the 
boundary  at  -z  has  the  same  group  of  particles  removed,  and  their  contribution  must  be 
added  to  (B9).  We  define  the  other  groups  of  boundary  particles  as 
S  '_  =  {particles  j,  such  that  Rzj-Rzi>L/ 2} 
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S~  =  {particles  j,  such  that  Rx  i  -  RX  J  >  L/ 2} 
Sx  =  {particles  j,  such  that  Rx  .  -  Rx  .  >  L  /  2  } 
S~  =  {particles  j,  such  that  Ry  .  -  Ry  .  >  LI2  } 


S+  =  {particles  j,  such  that  Ry  .-Ry  .>  LI  2} 

where  the  +  or  -  signifies  whether  the  contributions  from  the  particles  is  added  or 
subtracted  from  the  +boundary. 

Boundary  corrections  only  affect  and  C,'  , '!J  terms.  At  a  distance  far  from 

the  center  the  contribution  from  these  particles  on  the  boundary  can  be  derived 
analytically.  Far  from  the  center,  the  sum  over  cells  (cells  are  not  broken  into  subcells) 
can  be  written  as  an  integral  as  follows. 


Lcr"A,,,,®J+ 


sg+z  boundary 


sg-z  boundary 


L 


dAC, 


0,1  m 


K-@J 


+z  boundary 


-z  boundary 


(BIO) 


Using  the  off-center  expansion  of  Appendix  A  and  evaluating  the  above  integrals  yields 
the  correction  for  the  boundaries.  We  do  the  arithmetic  explicitly  for  the  +z  boundary  of 

the  evaluation  of  C^f  : 
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4 


+z  boundary 


N„L/2  N„L/2 


/  2  2 

-<C  f  dx  f  rfy  “y  2 - 7^,_m  arctan  V  *  ,arctan  - 

-ncl/2  -nc£/2  (N//2)  +x  +  y  1  I  NcZ/2  I 


(Bll) 


*^0»;  ,71.00  ^ 
L2  °’10V4^ 


J0ffl2^a; 


"  NCZ,  /  2  NcL/2 

j  dx  j  dy 

-N2./2  -tsLi/2 


a;NZ/2 

J  c 


(NcI/2) 


12  ,  2  ,  2 

1  +*  +y 


where  the  contributions  from  Cg  °°n’'oJs  and  are  taken  from  Appendix  A  and  the 

restriction  to  m  =  0  is  due  to  zero  contribution  from  the  integral  for  other  values  of  m. 

The  other  integrals  from  (BIO)  can  be  performed  similarly.  Contributions  from  the  x  and 
v  boundaries  yield  similar  expressions.  We  summarize  below  the  expressions  for  the 
boundary  corrections. 


=  -^o  when  j  e  5; 
+  c0  when  j  e  S. 


(B12a) 


K'°m  =  ±“/=  when  j  e  Sx 


+  ~j=  when  j  g  S~ 
~i—f=  when  j  e  .S’, 

V  2 


(B12b) 


+  i  —jL  when  j  e  .S\ 
v2 
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/e* = +ci when  j e  s: 

-  Cj  when  j  g 


/^0,1±1 


_£i_ 

2V2 


when  j  e  SA! 
when  j  e  S~ 
when  j  e  S' 
when  j  e  S~ 


n\,2±l,ij 

A'o.io 


=  ±^=  when  j  g  5; 


+  when  j  g  S; 


£ 

+  /  —jL  when  j  e  S 

a/2 


+ 

7 


£ 

-  *  when  j  e  S; 


K'J=Kl",J=-c2  when  j  e 

+  c2  when  j  g  Sz“ 

AuiY  =-|(±V3c,  -  2c2)when  j  g  SJ 
+  ^f±  V3c,  -  2c2  )  when  j  g  Sa. 
-7(^1  +  2c2)when  j  g  S,1, 


+  (V3cj  +  2c2  )  when  j  g  Sv 


(B12c) 


(B12d) 


(B12e) 


(B12f) 
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Kj±?’ij  =  ~(±  V3cj  +  2c2)when  j  e  St+ 

+  ^-(±  Scx  +  2c2  )when  j  e  S~ 

-  (V3c,  ±  2c2)when  j  e  S* 

+  ^-(V3q  ±  2c2)when  j  e  5“ 

with 

4;z»2 

C0  = - y- 

0  3Z2 

_  gj  [~  10V2;r(l  -  2cr)  5V6 
Cl_L2  3(4-5c)  +  (4  -  5cr) 

_  a)  5V2  lO^/l^l-cr) 

°2  ~~  Z?  (4  -  5cr)  (4 -5a) 

All  other  fi',  '"”  Jj  are  zero. 

Appendix  C.  Mechanical  properties  for  a  particle  pack  of  spherical  geometry 

The  mechanical  properties  of  a  spherical  pack  of  particles  can  be  estimated  from 
the  exact  solution  of  an  isolated  spherical  elastic  body  of  radius  a  with  shear  modulus  jue 
and  Poisson’s  ratio  ae  embedded  in  a  uniform  elastic  medium  with  shear  modulus  //  and 
Poisson’s  ratio  cr(Goodier,  1933).  For  uniaxial  tensile  displacement  with  strain  Ezz,  we 

can  compare  the  exact  Goodier  solution  to  the  multipole  solution  to  order  (r  is  the 

r 

distance  from  the  sphere  center)  to  relate  the  mechanical  properties  of  the  embedded 

sphere  to  the  multipole  moments.  In  spherical  coordinates  and  to  order  ,  the  Goodier 

r 

solution  for  displacement  in  the  region  far  from  the  sphere  is  given  by 


(B12g) 


(B13) 
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A  C( 5 -4a) 

(l  —  2cr)r 


u„  = — r+  /  x  — ^fcos(2#) 


ue  =-^sin(26>), 
r 


(Cl  a) 

(Clb) 


where 


A  = 


Bj\  +  a)o 


M~  Be 


(!  -  2o-g  X6  -  5^ )2//  +  (3  +  1 9ae  -  20 acre  )/ue 


(l  -  5cr)//  +  (8  -  1 0cr)/4 


(l-2o-j2//  +  (l  +  o-J/4 


(\-cr) 


+  2- 


1  +  cr, 

- --a 

1  +  a  1  j 


/4-(l-2<7> 


(l  -2a  )2//  +  (l  +  )//, 


c  _  E=  (l  +  5(l-2 crtfi- Me) 

4  (7  -  5cr)//  +  (8  -  1 0 cr)//g 

The  multipole  solution  for  the  pack  can  be  written  to  the  same  order  in  r  using  (3)  and  (9) 
for  the  displacement  field.  Because  the  B'0  hn  and  B[lm  are  zero  for  a  pack  in  equilibrium 


(i.e.,  forces  and  torques  on  each  particle  are  zero),  the  only  terms  that  contribute  to  — 

r~ 

are  the  B[  00  and  B[  lm  coefficients.  Only  the  m  =  0  tenn  contributes  in  B[  2m  for  a  large 

random  pack  since  the  m  ±  0  terms  imply  that  there  is  not  azimuthal  symmetry.  Keeping 
only  these  terms  it  is  a  simple  exercise  to  write  (3)  and  (9)  in  the  form  of  (Cl).  In  tenns 
of  the  multipole  coefficients  the  constants  A  and  C  are 


A  = 


C  = 


-  1  V  2p/  5(5 -4cr)  y  2p/ 

“  rV~2^aiBl,00  nr-(.  r\2^aiB  1,20 


-^4n  i  ‘  16V2^(4-5a)^T 

15(1  -2a) 


16V2tt(4  -5a) 


^jai  Bl,20 


(C2) 


(C3) 
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By  equating  A  and  C  in  (Cl),  (C2),  and  (C3),  one  can  solve  for  the  effective  mechanical 
properties  of  the  spherical  pack  of  particles  in  terms  of  the  multipole  coefficients  and  the 
mechanical  properties  of  the  binder  matrix. 


A,  = 


,_iZzMF 

(8-1  Oct) 

1 +  F 


A 


CT„  = 


G(2a  +  Ae )  +  (l  +  o~)a  ~  (l  ~  2crK 

g(4A  -  A  J  +  2(1  +  ct)//  +  (l  -  2 o)ne 


3  ZafBt 


20 


F  = 


E_2^i2tz[\  +  (7)c 


1,00 


G  = 


(C4a) 

(C4b) 


APPENDIX  D.  Mechanical  properties  for  an  infinite  pack 

In  this  appendix  we  derive  the  equations  that  allow  the  mechanical  properties  for 
an  infinite  pack  to  be  calculated  from  the  coefficients  of  the  multipole  expansion.  We  do 
not  limit  this  derivation  to  isotropic  packs.  In  general  the  effective  mechanical  properties 
of  a  medium  can  be  derived  from  the  average  stresses  and  strains  in  that  medium.  For  a 
material  made  up  of  rigid  particles  embedded  in  an  elastic  matrix,  the  average  stress  r.. 

(/,/  indices  range  over  x,  y,  and  z  values)  in  the  material  is  given  by  the  average  of  the 
sum  of  stresses  in  the  particles  and  the  binder  as 

<d» 

v  v  p 

where 
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(D2) 


Jy  =  J  TijdV  and  BU  =  y  I  Tydv  ■ 

Particle  p  Binder 

The  volume  of  the  system  under  consideration  is  V.  For  void-free  material  By  is  related 
to  the  average  strain  Ea  induced  in  the  medium  by 

B,J  =  +  K  +  Kh „  +  2  m(E°)  (D3) 

where  A  and  ju  are  the  Lame  constants  of  the  binder. 

To  evaluate  the  average  stress  in  the  material  we  must  separately  evaluate  the 
contributions  from  the  particles  and  binder  as  in  (Dl).  The  contribution  for  the  average 
stress  in  the  particles  is  derived  next. 

Equation  (D2)  for  the  volume  integral  of  r;/  in  particle  p  can  be  turned  into  a 
surface  integral  as  ( iz  component  shown): 


P 

I?  =  |  rrd  V  =  |  zi:dxdydz  =  J  dz  |  ri:dxdy 


Particle  p  Particle  p 


-an  Surfaces. 


p  p 

=  J  dz  J  r\ :  ■  x  •  da  =  J  dz  J  hi  •  x  •  da 


-ap  SurfaceS! 


-ap  Surfaces  2 


where  the  surface  5)  is  the  plane  surface  inside  the  sphere  perpendicular  to  the  z  axis  and 
passing  through  z,  S2  is  the  surface  of  the  sphere  for  all  points  greater  than  z,  and  ap  is  the 
radius  of  particle  p.  The  identity  on  the  right  hand  side  of  (D4)  is  derived  from  the 
equation  V  •  x  =  0  for  a  body  in  equilibrium  and  Gauss’s  theorem  (with  signs  adjusted  to 
the  orientation  of  da). 


0  =  |  itf  ■  x  ■  VdV  =  J  hi  -  T-da 

Volume  bounded  Surface 

by  Sj  +S2  Sx+S2 


Surface  Sj 


t  •  da  +  J  n,  ■  x  •  da 


(D5) 


Equation  (D4)  can  be  simplified  by  switching  the  order  of  the  z  and  surface  integrals,  and 
performing  the  z  integral  explicitly.  The  result  is 
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(D6a) 


Jt=  frizdV=  \zhrx-da 

Particle  p  Surface  of 

Particle  p 


where  z  is  the  distance  in  the  z  direction  from  the  center  of  the  particle  to  the  position  of 


the  surface  element  d  a  .  Similarly, 


IF  =  Jx  hf  ■  x  •  da  (D6b) 

Surface  of 
Particle  p 


and 


Ify  =  \y  nt  •  T  -da.  (D6c) 

Surface  of 
Particle  p 


It  is  now  a  straightforward,  if  tedious,  exercise  to  calculate  the  IF  in  terms  of  the 
coefficients  B[nm  of  the  multipole  expansion.  The  stress  is  calculated  from  the  general 


multipole  solution  for  the  displacement  field,  and  then  substituted  into  (D6).  The  results 
for  the  volume  integral  of  stresses  in  the  particles  are  as  follows. 

1  -  a 


Pz  =  pAyfx 


1  u  V  g1Bp 

1-2 V2(4-5 v)YP  X2\ 


(D7a) 


Ipxx  =  pAyfx 


1  -  a 
1-2(7 


Ya2Bp  +  5(1~(J)  y  2L  _  [j(Bp  Bp  ) 

ZmJuPL>\, oo  ^  nr/  \A^up\r^  V  2  v"i,22  ^  1,2-2 /. 


2^2(4 -5a) 


(D7b) 


lyy  =  PA-Jn 


1  -  a 
1  -  2  a 


I  alB‘M  +  „  4  , 1  a l  [s' „  +  Vf  (S' ,  +  S' _2  \ 


2V2(4-5  a) 


(!-a) 


Ipzy  =lpyz  +^2-i) 

^  =1%  -^2) 


(D7c) 


(D7d) 


(D7e) 


(D7f) 
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We  next  evaluate  the  average  stress  in  the  binder  using  (D3).  To  do  this  we  must 
calculate  the  average  strain  Ea  in  the  material.  The  average  zz  component  of  strain  is 
given  by  the  change  in  the  ±  z  wall  position  of  the  unit  cell  divided  by  its  length  as 


Kz  =  iuz  (*,  y,  Z  +  L)-  u,  (x,  y,  z))/L 

=£“+jlIsr>v/Cv 


t  frf , 

Lnm  p 


z  boundary 


Y, 


10 


77 inf  .  ^  X^  DP  oL'n'm 

^zz  “l~  r -  L' ,n'm  y 010 

Lj'sJQTT  L'n'm'  p 


I  z  boundary 


_  77 inf  V  DP  _2  J_ 

“  Bzz  +  jJ, 


572^(1 -2  q)|  5 


3(4-5a)  (4-5o-) 


P  2 

20ap 


(D8a) 


where  the  contributions  to  the  difference  of  the  displacement  over  a  periodic  cell  come 
only  from  the  far  field  strain  at  infinity  EmJ  and  the  boundary  corrections  (equations 


(B12)  and  (B13))  derived  in  Appendix  B.  Similarly,  it  can  be  shown  that  the  other 
components  of  average  strain  are  given  by 


J7a  =  pint  _ 

■*"'  XX  ^  XX  2)1^ 


1 


IX 


p  a2  - 

00  up 


2  If 


5^2ni}  -  2o-) 

3(4 -5a)  +  (4  -  5a) 


A  p 


(D8b) 


e:„  =e™  ~2 


yy  yy 


3  L 


3 


1,00 a  P 


2  u 


5^2tt{\  -  2a) 

3(4 -5a)  +  (4  -  5a) 


(D8c) 


Eazx  =  Eaxz  =  E™  +  \ 

AX  XA  XZ  -y  j 


5  5jf(l  -a) 

2V^(4-5a)  (4 -5a) 


(D8d) 


Fa  =  Fa  =  Finf  +  — 

nzy  nyz  nyz  +  jJ3 


5  5^(1- a) 

2V^(4-5a)  (4 -5a) 


+su-.) 


(D8e) 
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Y.KiKi-i-Kn)  (D8fl 

P 

Equations  A-6.5  and  A-6.6  can  be  substituted  into  equations  A-6.1  through  A-6.3  to 
calculate  the  average  stress  tensor  in  the  medium  from  the  coefficients  of  the  multipole 
expansion.  Using  the  average  strain  tensor  calculated  in  eq.  A-6.6  then  allows  the 
mechanical  properties  of  the  medium  to  be  extracted. 


Fa  =  Fa  =  Finf  +  _L 

^  xy  C xy  +  jj? 


5  5jf(l  -a) 

lyfniA -5ct)  (4-5cr) 


Appendix  E.  Translation  coefficients  of  Navier  multipole  solutions  valid  far  from 
the  surface  of  sphere  i 

In  Appendix  A  off-center  expansions  were  derived  for  the  multipole  solution  of 
sphere  j  on  the  surface  of  sphere  i  in  a  reference  frame  centered  at  i.  A  different  set  of 
expansions  is  needed  for  points  far  from  the  surface  of  sphere  i.  These  expansions  are 
valid  when  the  distance  r.  from  sphere  i  is  greater  than  the  separation  R  /7  between 

spheres  i  and  j.  The  derivation  is  very  similar  to  that  presented  Appendix  A.  For 
brevity,  we  present  only  the  results.  The  nomenclature  remains  the  same,  with  the 
exception  of  the  definition  of  A  which  is  defined  as 


?ji=n-Rji  ■ 


(El) 


The  off-center  multipole  expansion  for  the  harmonic  term  is 


(  ..  A 


L,+l 


\rji  J 


y,%  (fi, .  9,  >-££  Z  A  A'  <*, .  n )  A,  w .  9, ) 


(E2) 


rif  =0  L;  =72,  -1  722,  =-  72, 


where  the  expansion  coefficients  ( R  „ ,  ^  )  are  functions  only  of  the  relative  position 


Lj  ;nt  ,mi  V  ji  ?  / 


vector  R  u  and  the  distance  r.  from  sphere  i. 
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The  off-center  multipole  expansion  for  the  nonharmonic  term  is 


(2n.-iw«.+i)  f  -  f 

v  fl  y  v  fi  y 


1:  =0  L:  =flj  -1  Ifl:  =—tl: 


where  we  will  only  calculate  the  off-center  term  for  Lj  =  nj  -  1  since  the  nonhannonic 
tenns  are  zero  for  the  other  two  values  of  Lj. 

The  harmonic  term  coefficients  are  given  by  the  integral 


\d<p,\dj cosS,)  -r  («„?>,) 


and  the  nonharmonic  term  coefficients  by  the  integral 


-  !)>,(», +  1) 


Expressions  for  scalar  harmonic  function  off-center  expansions,  which  are  needed  in  (E4) 
and  (E5)  for  the  functions  of  (/y,  tf/;, r/y),  are  given  by  Varshalovich  (1988): 


^  t=  H2Z  +  1> 

ri+1  ^  \  (2  L)\ 


t^+ily  1}; 
(2L)\  h 


(2T)!  cij+lRj2L 


(2(A-L)+l)\  rt 


l _ JJ_ 

A+\ 


■  t{A  x  L  L]w 


Performing  the  integral  of  (E4)  for  the  harmonic  term  of  the  Navier  off-center  expansion, 
one  obtains 


-  aLj+lRL‘  Lj 

Lj  >nj  'mj  /  5  „  \  _  ,/  Lj’nj'mj  j  ji  y 

/;  5  uLl:ni,m,  L  + 1  ‘  / 


rmi  _l - y - F,  ,  (0  0)  ) 

m;  ^  Lt  +1  Lj  -Lj  ,rrij  -rrij  V  yt  ’  Jl  ' 


where  the  constant  coefficient  is  given  by 
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■ 

di!Z£  =  H)  V"'  V4^(2  Lj  +  1X2«,  +  lX2«y  + 1)  j 


(24)! 

(2(4  ~  Lj)  +  1))(2Z, .)! 


L, 


1  n, 


77 


\m,+\  -1  -m,-A 


4  4-4 


L, 


77 


m,  +  1  m,  -m,  -  m  ■  -  \ 


L, 


1  nj 


4  1  n, 


77 


\mi  0  -mijy...t  -j 
r  L  1  n  V 


4  4-4  4 

111 ,  III  ;  -  11  1 .  -  11  1 


\mi  ~  1  1  -'",7V 


4  4-4 


y  7VV 

4  1  ni 

i  J\mi  0  ~w7  7 

77 


m ,  + 1  - 1  -  111 


JJ 


4 


m,  - 1  m  j  -  mi  -  m .  + 1 


4  1  nj 

m  .  —  1  1  —  m , 


7 


j  vv-v  "vv 

The  nonharmonic  tenn’s  off-center  expansion  is  as  follows: 


=  -l)4Tf7l)  +  ^sr>(«„r,) 


2 

4 


+  e  J  J 

~  c'4  ,/M: 


4  V; 

4 


r  a .  ^ 

”7 

_ J_ 

_ 

v  4'  J 

4-1 


'^Ll-nJ-3,mJ-ml  i®  ji  ’  ^  /,  ) 


V  r,  / 


4  L,:n;,m; 


fR,  7ii+1 


V  ri  J 


pnfmi 
Lt  ;rij , m. 


4 


+  g4 


7  ’V 


4 


4 


V  ri  J 


Lf-n  :+\,m  j-rrii  V  ji  ’  ji 


(®«>®*)  • 


(E8) 


(E9) 


Detailed  evaluation  shows  that  the  coefficient  444  's  always  zero.  It  will  not  be 
written  out  explicitly. 

The  coefficient  is  given  by 


A,,. ,my  =r_n»,+mi+i  4  (2/7,  +  l)(2n7  - 1)(2»7  +  3) 

-(  )  (2nj-2)\ 

x  [  4  (4 > mj ;  4 > rh > 4 )  +  4  (4 , ;  4 > 4 , 4 ) + 73  (4 , 4 ;  4 > 4 > 4 )  ] 

where 


(E10) 
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f\ (n., m . ; Z( , nt , mj )  =  ■ 


L .  - 1  L.  -  n ,  -  2  n ,  + 1  ^ 

1  1  J  J 

.  /«,.  -  2  -  m,  +  in  + 1  -  m .  + 1  . 

V  4  1  J  J  / 


x  ^(Z;  +  -  2)(Z(  +  mi  -  1 )( Z;  +  mt  -  ~nj  -  in ;  -  2)(Z(  +  m,  -  /? ;  -  /« ;  - 1) 

+  ^(Z,.  -  /«,  +  1)(Z,  -  m,  -  iij  +  nij  - 1) 


■  2V(A  +  mi  ~  1XA  + '»,■  -  -  m,  -  i) 
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1  J  J 
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«  V 


V  1  V 


f - f  Z,.  - 1  Z,  -  nj  -  2  n  .  +  1  ^ 

V(A- -  -  2) 


.  ni  -  m ,  +  m ,  - 1  -  m ,  + 1  , 

V  1  1  J  J  j 


V(2A-2)! 
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v  m,  - 1  1  -  m, 


Aw/ -1  1  -m/y 
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z  - 1  z  -  /?.  -  2  /?.  +  r 

*  1  J  J 
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m ,  - 1  1  -  m , 


2V"V 
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(ElOa) 
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f3(nj,mj;Lt,ni,mi)  =  - 


Z  —  1  L-n-  2  n . +1  ^ 

1  1  J  J 
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The  coefficient  is  given  by 
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F2(nj,mj;Li,ni,mi)  =  -\ 
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pL,.  +  1)(2Z,  +  3)(2Z(.  -  2 nj  +  1)(2Z,.  -  2«.y  +  3) 


n;  + 1 


/«,  -  /«,  +  w ,  - 1  -  ntj  +  \j 


L:  +  1  Li-nj  nJ  +  1 N 

g2(nl,m,;Li,ni,mi)  =  ~{  \ 

ni  - 1  -  /h,  +  //z  ,.  + 1  -  m . 

v  1  1  J  J  / 

x  ^(Zf  -  /zz,  +  1)(Z(  -  m~  +  2)(Z,  +  m,.  - zz~  - m/.)(Li  +  mi  - n /  - m j  + 1) 

+  +  Y  + !)( A-  “  m,  ~  n,  + '«/  + 1) 

Z,  +1  Z.  —  zz n , .  + 1 


2a/(Z/  -  /«,  +  1)(Z,.  +  /«,  -  iij  -  mj  +  1) 


+  p,  +  Y  +  2  )(Z,.  -  m,  -  tij  +  m2) 
Z,  1  ti,  Y  zz ,  + 1  1  ii 


v  « 

Y..  + 1 


/zz,  -  /zz.  +  m ,  -  m 


1  J 


L.-n. 


v  y 

n  j  +  Y 


/zz,  +1  -  ni  +  m—  1  -  /zz 


v 


‘  j 


' i  J 


m ..  0  -  ni 


\  i 


y v  "v 


m ,  0  -  m 


JJ 


pLt  +  1)(2Z;  +  3)(2Z,.  -  2/i ,  +  1)(2Z.  -  2/z .  +  3) 
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'  Lj  +1  Lj  —  n .  n.  + 1  ^ 

gi(nj,mj;Li,ni,mi)  =  \ 

^  mi  -mi  +  mj+ 1  —  1  1 

x  V(A-  -  m,)(A  - m,  +  l)(Li  +  mi  ~ n,  ~  mj)(Li  +  mi  ~ nj  -  m,  +  0 

+  +  mi  +  2X4  - mi  -  n,  +  m,  + ') 


, - f  L.  +  1  L  -  n 

2 7(A  -  mi)(Li  +  »h  -  nj  -  ntj  + 1)1 


tij  + 1 


ymi  +  1  - mi  +  ntj  - m j  - 1^ 


■  ^L.  +  m,  +  3)(L  -  m,.  -  «y  +  /«,.) 


L  + 1 


L,  ~  « / 


nj  + 1 


.  /«,  +2  -  /«,  +  in - 1  —  m ,  —  1 , 

V  1  1  J  J  J 


L,  1  «, 

in,  +1  - 1  -  m 


V-i 


^  rij  + 1  1  rij .  ^ 

v  m  v  +  1  -1  -  m  7  y 


a 


pL,  +  1)(2Z,  +  3)(2Li  -  2nj  +  1)(2Z,  -  2«.y  +  3) 


(E12c) 
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